SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Processing basespace generated bam file in DESeq2 Rammohan Bioinformatics 1 02-05-2016 05:22 PM
TCGA RNAseq BAM Files fk566938 Bioinformatics 2 12-29-2015 09:27 AM
Calling SNPs from RNAseq bam files jmwhitha Bioinformatics 7 10-15-2014 05:16 AM
MarkDuplicates not working on CG BAM files? biscuit13161 Bioinformatics 0 04-29-2013 05:28 AM
"R Killed" when working with large BAM files mixter Bioinformatics 2 07-05-2010 12:47 AM

Reply
 
Thread Tools
Old 04-11-2016, 03:39 PM   #1
Rammohan
Member
 
Location: Toronto

Join Date: Dec 2015
Posts: 10
Default BAM files from RNAseq -Alignment (Basespace) not working with DESEq2 in R

Hello experts,

With my recent run in illumina Hiseq2500 I aligned the FASTAQ files generated in basespace to human reference genome (Homosapieans (PAR masked)/hg19 (Refseq) using the RNA-seq Alignment App.I used Tophat (Bowtie2) aligner for the same. except for "NO" to the stranded information, No other parameters were selected.
I downloaded the BAM files generated and used them in the Bioconductor RNAseq workflow (Using DESeq2)
http://www.bioconductor.org/help/workflows/rnaseqGene/

to create a summarized object I used "Homo_sapiens.GRCh37.75.gtf"

below is the script:

(txdb <- makeTxDbFromGFF("Homo_sapiens.GRCh37.75.gtf", format="gtf", circ_seqs=character()))
(ebg <- exonsBy(txdb, by="gene"))

se <- summarizeOverlaps(features=ebg,reads=bamfiles,mode="Union",singleEnd=FALSE,ignore.strand=TRUE,fragments=TRUE )

I am getting the following error message and am unable to get any count

Warning message:
In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)

Any suggestions?

Ram
Rammohan is offline   Reply With Quote
Old 04-12-2016, 01:34 AM   #2
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 122
Default

Hi Ram,

could you check if the chromosomes' name in the bam files are the same as the ones in the gtf?

Code:
samtools view -H my_bam.bam | head
and
Code:
head Homo_sapiens.GRCh37.75.gtf
.

Often, one annotation is 'chr1', another 'Chr1', and sometimes just '1'. If you are mixing these up, no program or tool will find an overlap.

Cheers,

Michael
Michael.Ante is offline   Reply With Quote
Old 04-12-2016, 07:00 AM   #3
Rammohan
Member
 
Location: Toronto

Join Date: Dec 2015
Posts: 10
Default

Thanks a lot Michael,

below is the output of the codes you suggested:

[rshukla@dev01 test]$ samtools view -H Ctrl13058.bam | head
@HD VN:1.0 SO:coordinate
@RG ID:0 SM:Ctrl13058-Trimmed PI:50
@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663

[rshukla@dev01 test]$ head Homo_sapiens.GRCh37.75.gtf
#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
1 pseudogene gene 11869 14412 . + . gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene";
1 processed_transcript transcript 11869 14409 . + . gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; g ene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana";
1 processed_transcript exon 11869 12227 . + . gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; exon_numb er "1"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; exon _id "ENSE00002234944";
1 processed_transcript exon 12613 12721 . + . gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; exon_numb er "2"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; exon _id "ENSE00003582793";
1 processed_transcript exon 13221 14409 . + . gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; exon_numb er "3"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; exon _id "ENSE00002312635";


The format in .bam file is chr1 where as in .gtf I guess it is just the number.

Do you have any suggestion how to proceed further, I am very new to all these.

Thanks again,

Ram
Rammohan is offline   Reply With Quote
Old 04-12-2016, 07:09 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,992
Default

Can you download a corresponding GTF file from BaseSpace? That should sort this problem out.

It should be possible to change the GTF file you have (from 1 to chr1 etc) but if the above works then you would know that you have a file that exactly corresponds to the genome used for alignment. Using an incorrect GTF file will lead to erroneous results.

Edit: If changing 1 to chr 1 is the preferred option then use the solution suggested by @vivek in this thread.

Copying it here for reference (use your file names) :

Code:
$ awk '{if($0 !~ /^#/) print "chr"$0; else print $0}' no_chr.file> with_chr.file

Last edited by GenoMax; 04-12-2016 at 07:22 AM.
GenoMax is offline   Reply With Quote
Old 04-12-2016, 09:24 AM   #5
Rammohan
Member
 
Location: Toronto

Join Date: Dec 2015
Posts: 10
Default

Thanks GenoMax and Micheal

I tried the other reference genome with same chromosome annotation but it still did not work

Below is the error message

> se <- summarizeOverlaps(features=ebg, reads=bamfiles,mode="Union",singleEnd=FALSE,ignore.strand=TRUE,fragments=TRUE,suppressWarnings())
Error in validObject(.Object) :
invalid class “SummarizedExperiment0” object: 'assays' nrow differs from 'mcols' nrow


the heads of both the bam and the gtf files are as below:

[rshukla@dev01 test]$ samtools view -H Ctrl13058.bam | head
@HD VN:1.0 SO:coordinate
@RG ID:0 SM:Ctrl13058-Trimmed PI:50
@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
[rshukla@dev01 test]$ head Human_hg19.gtf
chr1 hg19_knownGene exon 11874 12227 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3";
chr1 hg19_knownGene exon 12613 12721 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3";
chr1 hg19_knownGene exon 13221 14409 0.000000 + . gene_id "uc001aaa.3"; transcript_id "uc001aaa.3";
chr1 hg19_knownGene exon 11874 12227 0.000000 + . gene_id "uc010nxr.1"; transcript_id "uc010nxr.1";
chr1 hg19_knownGene exon 12646 12697 0.000000 + . gene_id "uc010nxr.1"; transcript_id "uc010nxr.1";
chr1 hg19_knownGene exon 13221 14409 0.000000 + . gene_id "uc010nxr.1"; transcript_id "uc010nxr.1";
chr1 hg19_knownGene start_codon 12190 12192 0.000000 + . gene_id "uc010nxq.1"; transcript_id "uc010nxq.1";
chr1 hg19_knownGene CDS 12190 12227 0.000000 + 0 gene_id "uc010nxq.1"; transcript_id "uc010nxq.1";
chr1 hg19_knownGene exon 11874 12227 0.000000 + . gene_id "uc010nxq.1"; transcript_id "uc010nxq.1";
chr1 hg19_knownGene CDS 12595 12721 0.000000 + 1 gene_id "uc010nxq.1"; transcript_id "uc010nxq.1";
[rshukla@dev01 test]$


As per Genomaxs suggestion I am also downloading reference genome from basespace, surprisingly it is 42.3GB !

any further suggestions?

Ram
Rammohan is offline   Reply With Quote
Old 04-12-2016, 10:41 AM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,992
Default

@Ram you should not need to get the reference genome. The error above is not related to the original issue but an error with your DESeq2 command.
GenoMax is offline   Reply With Quote
Old 04-12-2016, 12:11 PM   #7
Rammohan
Member
 
Location: Toronto

Join Date: Dec 2015
Posts: 10
Default

Thanks, what changes should I make to the command?
Rammohan is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:58 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO