This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract specific fragments from SAM and count them

Hi, friends. I have a .sam file & .vcf SNPs file. I want to extract the sequnece fragments which overlaping given heterozygous SNPs sites. Then I want to classify the fragments by the different base on SNPs site and count the number N0 & NX. Can anybody tell me how to code the shell script or which software can do that.

snp sam
• 4,222 views • link
updated by BioinfGuru BioinfGuru • written by Jeason Rad • Jeason Rad

2 answers

If you want to extract reads at a given position you can use samtools (ex chromosome 5 position 513931) :

samtools view -b file.bam chr_5:513931-513932 > reads_position.bam

with a list of positions extracted from your vcf, you could loop over it to extract each position :

for i in $(cat list_positions_SNP); do samtools view -b file.bam ${i} > reads_${i}.bam;done
• 0 views • link
written by guillaume.rbt guillaume.rbt

Thx, friend. But I puzzle more about how to number No and Nx basing on the different base at SNVs site. I had uploaded a sketch map about it. What I want is not the total fragments on given region, but the number of No & Nx.

• 0 views • link
written by Jeason Rad • Jeason Rad

maybe this information is already in your vcf file, you should look for the field "AD" (allelic depth), which should report the coverage of both alt and ref alleles

• 0 views • link
written by guillaume.rbt guillaume.rbt

yes, AD is what I want. But my ultimate aim is to classify alt alleles fragments by their start or end sequencing sites (so as ref). enter image description here I mean I need to keep classifying and number the fragments by their positions (if fragment start or end in binding site). In the end, figure out No, Nx, no, nx; Any suggestions will be appreciated.

• 0 views • link
written by Jeason Rad • Jeason Rad

I use this command-line getting each fragment's pos/start/end/sequence in BAM.

bedtools bamtobed -i reads.bam -bedpe | awk -v OFS="\t" '{if($9=="+"){print $1,$2+4,$6+4,$10}else if($9=="-"){print $1,$2-5,$6-5,$10}}' > fragments.bed
• 0 views • link
written by Jeason Rad • Jeason Rad

Ok I understand. Indeed it's not that trivial. You have to parse the CIGAR string to know if your read is alt or ref at this position.

This topic has already been discussed here : Extracting Reads Containing A Specific Variant From A Bam File

It seems that the person that have posted the question have written a script that does the trick, maybe you can ask him.

• 0 views • link
written by guillaume.rbt guillaume.rbt

Convert to position sorted bed files and use bedtools intersect

• 0 views • link
written by BioinfGuru BioinfGuru

But if I turn SAM to BED, I don't know what the base on SNPs site. I had uploaded a sketch map about it.

• 0 views • link
written by Jeason Rad • Jeason Rad

could you re-write "I don't know what the base on SNPs site" please? It makes no sense.

• 0 views • link
written by BioinfGuru BioinfGuru

Log in to answer this question.

More posts like this