Question: Error in raw count file of .sam after alignment on GRCh38
0
gravatar for Mirza.Jawad
20 months ago by
Virtual University of Pakistan
Mirza.Jawad0 wrote:

Hello! I downloaded the indexes of GRCh38 human genome from ftp://ftp.ccb.jhu.edu/pub/data/bowtie_indexes/ and aligned my pair ended reads on it. after getting .bam, i converted it into .sam by using

**"samtools view -h -o output.sam input.bam"**

then i run

**"htseq-count output.sam Homo_sapiens.GRCh37.75.gtf "**

then I sorted the .sam file by using command

**"samtools sort -O sam -T 13.sam -o 13.sort.sam 13.sam"**

but it results like

...........................................................................................................................................................
Warning: Read SRR993713.6118094 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
Warning: Read SRR993713.122141 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
Warning: Read SRR993713.10023628 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
...........................................................................................................................................................

and then it shows the result with ensemble gene ids with 0 counts. please tell me whats wrong with it. Am i using the wrong gtf file. or the way of sorting my file is wrong?

ADD COMMENTlink modified 20 months ago by Macspider2.6k • written 20 months ago by Mirza.Jawad0
2

.bam, i converted it into .sam

why ?

ADD REPLYlink written 20 months ago by Pierre Lindenbaum115k

I learnt the pipeline from a person he said that u need to align then i have to get .sam formate but tophat gave me result in bam formate then i thought i should convert it into sam formate. what is the desired output formate of tophat tp get raw counts?

ADD REPLYlink written 20 months ago by Mirza.Jawad0

You can use featureCounts to get a count matrix directly from your BAM files.

ADD REPLYlink written 20 months ago by genomax59k

Its not working at all :(

ADD REPLYlink written 20 months ago by Mirza.Jawad0

It is not working. I spent 3 hours on it :(

ADD REPLYlink written 20 months ago by Mirza.Jawad0

What is not working?

ADD REPLYlink written 20 months ago by genomax59k
2

Am i using the wrong gtf file.

Yes you are. As @Ram pointed out, why are you mixing genome builds. This is a surefire way of getting erroneous results. Consider yourself fortunate that you got 0 counts, otherwise you may have continued on your merry way with completely incorrect data.

ADD REPLYlink written 20 months ago by genomax59k

then whats the solution ?

ADD REPLYlink written 20 months ago by Mirza.Jawad0
1

Use a GTF file for GRCh38 (since you used indexes for that build). Did you use Ensembl? The GTF files for that build can be found here.

ADD REPLYlink modified 20 months ago • written 20 months ago by genomax59k
1

Why download the GRCh38 indexes and use them with the GRCh37 gtf file? Am I missing something here?

ADD REPLYlink written 20 months ago by RamRS19k
3
gravatar for Macspider
20 months ago by
Macspider2.6k
Vienna - BOKU
Macspider2.6k wrote:

When htseq-count shows that error, it means that some reads which were supposed to be paired end didn't have a mate. When you run htseq-count with paired end mode, it wants all the reads to have a mate, so you have to filter keeping only the proper pairs first (correct orientation and within insert size).

samtools view -b -h -f 0x2 input.bam > output.bam

This will account for those warnings about reads without mates.

If your final counts file doesn't have any counts inside, you could be using the wrong version of the GTF file (you have to use the one which has been generated on the assembly that you mapped the reads against). As @genomax2 said, "Use a GTF file for GRCh38 (since you used indexes for that build)".

ADD COMMENTlink written 20 months ago by Macspider2.6k

I am getting sam result after using samtools view -b -h -f 0x2 input.bam > output.bam and after using htseq-count for Homo_sapiens.GRCh38.88.chr.gtf

ADD REPLYlink written 20 months ago by Mirza.Jawad0

An important indicator of how good your results are is in the lines which have an underscore before ("_"), in the htseq-count output. Try grepping them, you'll see the reads that haven't been assigned. If that number is extremely high, there are still some problems.

ADD REPLYlink written 20 months ago by Macspider2.6k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1181 users visited in the last hour