samtools complex question

Question

[biotech]

I need to know if there is any read mapping to certain region in an alignment file (.sam). In fact, I have a list of 20 regions I'm interested. I'm talking about the third column in the sam file. Also, the reads I want to filter must map with an alignment score greater than a certain value.

I'm starting from 6 sam files.

Maybe a custom script is better than samtools?

Answer 1

[Pierre Lindenbaum]

if there is any read mapping to certain region in an alignment file (.sam)

Samtools And Region List

Also, the reads I want to filter must map with an alignment score greater than a certain value.

Filtering A Sam File For Quality Scores

[biotech]

Thanks Pierre. I just need all reads mapping to certain chromosome. Would this ease the process? Not looking for coordinates.

[biotech]

Am I right with this shortcut? I have no output, shell running... Possibly due to large *sam files.

for file in ./*sam; do echo $file; samtools view -S -F 4 -q 30 $file; done | grep -x -f ./GBK/GBK_homologues/M6_f0_0taxa_algCOG_e0_/specifics.locus

[John]

Note that this would filter out unmapped reads but not mapped mates, giving a file with singletons and pairs.

[biotech]

And this one? Still no output

awk -F"\t" '{if ($5 > 30) {print $0}}' *.sam | grep -x -f ./GBK/GBK_homologues/M6_f0_0taxa_algCOG_e0_/specifics.locus

[Devon Ryan]

Please actually read the first link that Pierre posted. Grep should never be needed here.

[biotech]

Is this more or less what you suggested? 'specifics.locus' contains my interesting bed file first columns IDs that I'm trying to filter.

for file in *sam.bam.sorted.bam
do
    bamToBed -i $file | awk -F"\t" '{if ($5 > 30) {print $0}}' > $file.bed
    grep -x -f ./GBK/GBK_homologues/M6_f0_0taxa_algCOG_e0_/specifics.locus $file.bed &
done

[Devon Ryan]

BTW, this would be much faster if you made a sorted BAM file.