1
$\begingroup$

I want to compute the depth of coverage only for specific intervals in phase 3, 1000 genomes project.I have not worked with 1000 genomes project before, so a bit unfamiliar with it. I do not want to download all of the bam files for the entire genomic region, just those intervals which I am interested in them. I guess I can use this (I found from biostars forum)

samtools -bu view 'http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam' "2:1000000-2000000" | (...)

However, I am just a bit confused. First of all, when I try the code above directly from the command line, I get an error message "No such file or directory"! Am I doing something wrong? or should I provide a specific path? Second, I want to run it for the mapped files for all individuals: *.bam”. through a loop.

Can anyone help me to solve my problem?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

The command you're looking for is:

samtools view -c http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam 2:1000000-2000000

For a series of intervals, you'd be better off scripting something with pysam:

import pysam
bam = pysam.AlignmentFile("some file.bam")
regions = open("your regions.bed")
for region in regions:
    chrom, start, end = region.split("\t")
     bam.count(chrom, int(start), int(end))  # Do something with the result

You might see if there are bigWig files available, they'd be faster to process.

Improve this answer
$\endgroup$
5
  • $\begingroup$ ,thanks. I could rum the sam tools command for a single file. Thanks for the pysam, I will try it. But do you know how can I run a samtools command through a loop? just to remind that I do not have bam files downloaded. $\endgroup$
    – Anna1364
    Apr 23, 2018 at 21:18
  • $\begingroup$ You can do for loops in bash. It can be pretty painful to script anything complicated in bash, though. $\endgroup$
    – Devon Ryan
    Apr 23, 2018 at 21:20
  • $\begingroup$ yes, I know the for loop solution. my problem is that the bam files for each individual are in separate folders. I am trying this: for file in ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data*/alignment/*.bam; do /data/programs/samtools-1.3.1/samtools view -c "${file}" 2:1000000-2000000 done but get an error message. $\endgroup$
    – Anna1364
    Apr 23, 2018 at 21:46
  • $\begingroup$ You can't wild-card URLs, you'll need to provide them as a list. $\endgroup$
    – Devon Ryan
    Apr 23, 2018 at 22:17
  • $\begingroup$ @Anna1364 please add all this information to your question. We can't answer you if you don't tell us all the details. Give us a list of the files you want to read and explain what you need to do with those files. $\endgroup$
    – terdon ♦
    Apr 24, 2018 at 8:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.