Question: Extracting reads from given region
0
gravatar for KVC_bioinfo
23 months ago by
KVC_bioinfo410
Boston
KVC_bioinfo410 wrote:

Hello all,

I am extracting the reads from bam file which are falling in particular region using the following command.

samtools view -h input.bam "Chr10:18000-45500" > output.bam

However, now I want to extract the reads in sam manner but my reference is human genome primary assembly where the name of fasta reads starts like:

NC_000001.10 Homo sapiens chromosome 1, GRCh37.p13 Primary Assembly

How do I now modify command to do the sam task?

Thank you in advance.

bam • 1.0k views
ADD COMMENTlink modified 23 months ago by Pierre Lindenbaum124k • written 23 months ago by KVC_bioinfo410
3
gravatar for Pierre Lindenbaum
23 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum124k wrote:

the fasta should be indexed ignoring everything after the first whitespace. So your command should be something like:

samtools view -h input.bam "NC_000001.10:123-456" > output.sam
ADD COMMENTlink modified 23 months ago • written 23 months ago by Pierre Lindenbaum124k

Thank you. However, I want to extract reads from Chr 10 with the given coordinates above.

ADD REPLYlink written 23 months ago by KVC_bioinfo410
1

try looking for some thing like: NC_000010 in header:) $ grep ">" <input.fasta> (if your fasta is gizpped, use zgrep instead of grep). You would get all the headers in reference file. Then look for an entry that starts with NC_000010. Use that to extract the information. Command would look like:

samtools view -h input.bam "NC_000010.xx:123-456" > output.sam
ADD REPLYlink modified 23 months ago • written 23 months ago by cpad011212k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1742 users visited in the last hour