Question: Extracting reads from given region
0
gravatar for KVC_bioinfo
2.8 years ago by
KVC_bioinfo510
Boston
KVC_bioinfo510 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.4k views
ADD COMMENTlink modified 2.8 years ago by Pierre Lindenbaum130k • written 2.8 years ago by KVC_bioinfo510
3
gravatar for Pierre Lindenbaum
2.8 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum130k 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 2.8 years ago • written 2.8 years ago by Pierre Lindenbaum130k

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

ADD REPLYlink written 2.8 years ago by KVC_bioinfo510
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 2.8 years ago • written 2.8 years ago by cpad011214k
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: 1690 users visited in the last hour