This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Sequences information from BAM file

Hello, I'm working with paired end rna-seq data so two fastq files produced 1 BAM file. I know col 10 represent the seq and col 11 have its ascii code for quality. If the original seq length 100bp when I extracted some seq from BAM file I noticed some reads of length LESS than the original length. For example in one location one seq has cigar (40M, which means 40 bps matched) my question is: what happened to the rest of the seq (60bps)?? Another point, in this case we have paired end data which means two seq one from each end but BAM file has only one. This one included in BAM file to which end belongs??

Many thanks..

rna-seq next-gen alignment
• 3,631 views • link
updated by siqizhao2018 • siqizhao2018 • written by Lolla • Lolla

Did you map locally or end-to-end? Which aligner did you use?

• 0 views • link
written by Matteo Schiavinato Matteo Schiavinato

end-to-end using STAR

• 0 views • link
written by Lolla • Lolla

Are you 100% sure? The CIGAR you reported is compatible with a local mapping. What is your command, can you post it here?

For what concerns the absence of the paired mate from the BAM: if you used an option to leave out unmapped scores, perhaps that's the reason.

To know which end it does correspond to, you have to see the bitwise flag. Reads on the reverse strand will have the 16 bit.

What is the bitwise flag for the paired read that doesn't have the mate in the file? Did you check it on https://broadinstitute.github.io/picard/explain-flags.html ?

• 0 views • link
written by Matteo Schiavinato Matteo Schiavinato

Thanks Macspider; I confirmed trimmed happened prior the alignment, so this is solved now..

• 0 views • link
written by Lolla • Lolla

Then you can mark the thread as closed, so it doesn't show up in the "open threads" section. ;) Good luck!

• 0 views • link
written by Matteo Schiavinato Matteo Schiavinato

Thank you, I'm still confused about some fields in BAM file..

• 0 views • link
written by Lolla • Lolla

4 answers

it could be some HARD clipping in your bam ... Difference between Hard Clip and Soft Clip in Samtools CIGAR string or you reads could have been trimmed by a tool like cutadapt

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

There are basically two reasons that your reads would be trimmed like that...for quality reasons, or because the read read all the way through to adapter. If it's adapter, you definitely want that trimmed away. If it's for quality reasons..you probably want it trimmed away. You'll have to ask whoever did the alignment if they did a quality trim before aligning.

• 0 views • link
written by swbarnes2 swbarnes2

Another point, in this case we have paired end data which means two seq one from each end but BAM file has only one. This one included in BAM file to which end belongs??

Some aligners trim read headers after first whitespace when writing them to alignments, which removes Illumina identifiers signifying read number. Sounds like that applies in your case. So there is no easy way to figure out what end of the fragment the read came from.

• 0 views • link
written by GenoMax GenoMax

hi,I want to split paired end read as length distribution like (0-50,50-150, >150),How can I do it? when I using samtools extract reads(samtools view -h Epi2_input_rmdup.bam | awk 'length($10) > 100 && length($10) < 200|| $1 ~ /^@/' | samtools view -bS - > 100_200qc_rm/50_Epi2_input_rmdup.bam),I found all reads length no more than 100bp, but when I draw paired reads fragment size stastics using deeptools bamPEFragmentSize,it show length 0-500bp.

• 1 views • link
written by siqizhao2018 • siqizhao2018

Log in to answer this question.

More posts like this