samtools specific region

[dnusol]

I have downloaded all fasta files for human chromosomes and merged them with cat. In order to use the "region" option in samtools view command, do I have to specify the exact header of the fasta file?

so for example if I want to analyse only a specific region in chromosome 12:

$ samtools view gi|89161190|ref|NC_000012.10|NC_000012 Homo sapiens chromosome 12, reference assembly, complete sequence:1000000-2000000

I guess it's better then to modify the header to something like ">chr12".

Any suggestions?

Thanks

[krobison]

I would modify the header; if you don't you'll always need to put quotes around & type that whole mess. I don't believe in any case you would need the description information -- just the text up to the first space (but my samtools-indexed FASTA files lack descriptions, so I can't be sure)

How did you convert FASTA to SAM? Or did you mean to say "samtools index"?

[bioinfosm]

definitely, some sed commands to remove the ugly string and put >chr10 as header. Helps later on as well, when using UCSC or IGV!

[dnusol]

Thanks for your help

I indexed the fasta reference using bwa, but before doing that I changed the headers using sed. If after aligning against the whole genome, just wanted to get subalignments for a specific region, I have to use the samtools view command with the region option, right?

[dnusol]

Actually, just found out that the samtools view command does not work with the "region" option unless you feed an indexed BAM file, or so it seems:

$ samtools view -uS /s_1/s_1.sam.gz chr6:136000000:146000000 | ./samtools sort - /s_1/s_1
[samopen] SAM header is present: 25 sequences.
[main_samview] random alignment retrieval only works for indexed BAM files.

any suggestions?

[dnusol]

I will answer myself.

I first sorted and then indexed the BAM file. Then the region option seems to work.

D.

[krobison]

You just need to index (with "samtools index") the bam file that "samtools sort" generated -- then you are off to the races

[bioinfosm]

I think its sort, and then index, but perhaps it works in the other order as well

[genbio64]

Does anyone know of a way to use Samtools to split off individual chromosomes for easier viewing? If I attempt to use samtools view <myfile.bam> chr1 I receive a bam file that I can no longer view because it does not associate with the indexed file it was derived from.

[nilshomer]

Does anyone know of a way to use Samtools to split off individual chromosomes for easier viewing? If I attempt to use samtools view <myfile.bam> chr1 I receive a bam file that I can no longer view because it does not associate with the indexed file it was derived from.

Could you post the error message when it complains?
Try:

rm myfile.bam.bai
samtools index myfile.bam
samtools view -b myfile.bam chr1 > myfile.chr1.bam
samtools view -b myfile.bam chr2 > myfile.chr2.bam
...

[genbio64]

@nilshomer
I'll try that.

thanks