This is a test version of Biostars. For the public version, visit https://www.biostars.org.
filtering alignned seqeunces

Hi, I have used ONT minimap to align combined reference sequences (fasta file combined together in one file but separated by header " >" for each reference sequence). Then, I got one aligned sam file . My question is I dont know how to filter out my aligned data to each of the reference sequence. Details are below: I have 2 reference sequences in fasta format (lets call it A and B) And one aligned data in sam format. I want to see how many sequences align to each referenec DNA sequence (A and B)

Do I use the awk comand

alignment
• 1,325 views • link
updated by colindaven colindaven • written by priscilladraj • priscilladraj

I tried this 

awk '{if($3 == "A.fasta") { print $0 }}' combineA_aln.sam > A_aln_reads.file

wc -l A_aln_reads.file

to pull out A

awk '{if($3 == "B.fasta") { print $0 }}' combineB_aln.sam > B_aln_reads.file

wc -l B_aln_reads.file

to pull out B

but its not working. Can I know what is wrong. I am really new to comand lines in Linux

• 0 views • link
updated by lieven.sterck lieven.sterck • written by priscilladraj • priscilladraj

can you add what is not working or why you think it's not working?

• 0 views • link
written by lieven.sterck lieven.sterck

1 answer

A more standard way to do this would be

  1. convert SAM to BAM with BAI index
  2. run samtools idxstats x.bam
  3. alternative: run samtools stats x.bam

This would give you the stats including number of reads aligned to each reference as a text file. Very useful if you have many references.

• 0 views • link
written by colindaven colindaven

Log in to answer this question.

More posts like this