Question: Raw Read count from bam file
0
gravatar for Bioinfonext
3.3 years ago by
Bioinfonext200
Korea
Bioinfonext200 wrote:

I mapped pair end reads to CDS transcripts using BWA.

If I do sorting of bam fine by not using -n flag than its worting:

For sorting: /home/yog/software/samtools-1.3.1/samtools sort 216_5W_Ca1.bam -O BAM -o 216_5W_Ca1.sort.bam

For indexing: /home/yog/software/samtools-1.3.1/samtools index 216_5W_Ca1.sort.bam

For read count: /home/yog/software/samtools-1.3.1/samtools idxstats 216_5W_Ca1.sort.bam > readcount[bam_idxstats]

But read count I am getting like this:

Rs025080 1341 239 27

Rs035250 621 0 0

Rs035280 408 0 0

Rs035290 318 0 0

Rs035300 456 87 0

Why there is three count column for each transcripts.

rna-seq • 3.2k views
ADD COMMENTlink modified 3.3 years ago by Jorge Amigo11k • written 3.3 years ago by Bioinfonext200
1

I quick look into the manual gives you: "Retrieve and print stats in the index file. The output is TAB-delimited with each line consisting of reference sequence name, sequence length, # mapped reads and # unmapped reads".

ADD REPLYlink written 3.3 years ago by ATpoint32k

How do you have # unmapped reads per chromosome xD

ADD REPLYlink written 3.3 years ago by John12k
1

Here are two reasons: Samtools Idxstats

ADD REPLYlink written 3.3 years ago by ATpoint32k

Yes the 'Read is unmapped but has the chromosome of it's mapped pair" was my first thought, but then the sum of mapped/unmapped per chromosome would be an even number - but for the last entry it is not, so i don't think that's what's going on.

I didn't know about the BWA concatenating contigs though, that's new to me - thanks for the info :) And that could certainly be what's causing it. I also didn't know the SAM spec calls for unmapped pairs to be on the same chromosome - I always thought they came at the end with all the other unmapped reads. Two biologically meaningless but bioinformatically important bits of information today :)

ADD REPLYlink written 3.3 years ago by John12k

It stands for "number" in many manuals.

ADD REPLYlink written 3.3 years ago by Macspider3.0k
2

That's not what John meant, but rather "if a read is unmapped - how do you know it belongs to a certain chromosome" ;-)

ADD REPLYlink written 3.3 years ago by WouterDeCoster43k

BWA for instance assigns both the unmapped flag and a position. In case the read maps somewhere but is not fulfilling the criteria of being reported as mapped (e.g. too many mismatches).

ADD REPLYlink modified 3.3 years ago • written 3.3 years ago by michael.ante3.6k

So the read "almost" maps there :p

ADD REPLYlink written 3.3 years ago by WouterDeCoster43k

Here I map raw reads to a CDS transcripts of a draft genome. So there are around 30% reads of a paired-end library do not map to CDS.

But I am not sure how number were given of unmapped reads according to CDS transcripts ID.

ADD REPLYlink written 3.3 years ago by Bioinfonext200

Could you reformulate this question in a more explicit fashion? I am not sure if I got the point. I also have a question: when you say:

I map raw reads to a CDS transcripts

Do you mean against a transcriptome, or against a genome with a GFF/GTF/GFF3 file that indicates the CDS of the transcripts, and you are considering only those?

ADD REPLYlink written 3.3 years ago by Macspider3.0k

I took CDS from this database. but its genome is not completely sequenced. it's a draft genome.

http://radish-genome.org/Data_resource/

ADD REPLYlink written 3.3 years ago by Bioinfonext200

As far as I know, the notation CDS defines those regions of the exons that are translated, because exons also have UTRs. So you might have some RNASeq reads that belong to UTRs and are therefore not present in your CDS sequences. Could this be the case?

ADD REPLYlink modified 3.3 years ago • written 3.3 years ago by Macspider3.0k
0
gravatar for Vitis
3.3 years ago by
Vitis2.3k
New York
Vitis2.3k wrote:

I think to get read mapping stats you should use "samtools flagstat" or for the newer version of samtools, "samtools stats".

ADD COMMENTlink written 3.3 years ago by Vitis2.3k
0
gravatar for Jorge Amigo
3.3 years ago by
Jorge Amigo11k
Santiago de Compostela, Spain
Jorge Amigo11k wrote:

to get read counts you can use -c option, plus all the other options you may want to use to filter the reads to count, like -F4 to count mapped reads only:

samtools view -c input.bam
ADD COMMENTlink written 3.3 years ago by Jorge Amigo11k
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: 2227 users visited in the last hour