extracting identical reads from bam file

[Schelarina]

Hello,
I have an alignment file of smallRNAs in bam format.
How can I extact from this file all the identical reads in a single sequence but mantaining the reads count?
Thank you for the help!

[Richard Finney]

Are you trying to match a specific sequence ?

samtools view x.bam | grep TTTTCTGCCTGTTGGGCTGGAG | awk '{print $10}' | uniq -c


All reads with same sequence as another read, try this ...

samtools view x.bam | awk '{print $10}' | sort --buffer-size=20G | uniq -c | awk '{if ($1!=1) print $0}'


where x.bam is you bam file.
fine tune the --buffer-size parameter to sort

[swbarnes2]

It would take longer, but swap the second awk in Richard's command line with

sort -nr

To get the most common reads in order, starting with the most abundant. Tack a

| head -n 100

to get only the top 100

[Schelarina]

thank you! it worked!
can I do the same from a fastq file?

[Richard Finney]

Yes, you might want to use awk to only print every 4th line.
Note 1) the "mod" operator and 2) "line count" intrinsic variable in awk.
Perl/python/C/java if you prefer can address the issue of filtering for only the sequence also.

[Schelarina]

After I extracted the identical reads in a single sequence from the bam file I aligned them again to the genome. When I use igv to visualize the alignment now all the sequences are mapping in sense orientation.. even those sequences that are supposed to be antisense to the genome are shown in sense orientation. Why is that?

[Richard Finney]

The sequence for the read in a bam files may be reverse complemented to align to the reference.
Reads are supposed to be properly noted as reversed in the bitwise flags field in a line/entry of sam/bam.

You may wish to interrogate this flag for special processing.