Hot answers tagged

23

You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but since no one has added NCBI names (yet, maybe I'll add them now) you're currently out of luck there. Anyway, the quickest way to do what you want is to samtools ...


18

It's not really possible to convert bam to vcf. bam is a mapping file, it does not contain the information about variants, this information needs to be inferred in process called variant calling. I find important to mention that it's not just a different format of the same thing. How to call variants (a vcf file) from mapped reads (a bam file) is very broad ...


17

Whenever you want to save space (this can be a substantial savings). Until quite recently (samtools/htslib 1.7), only CRAM supported long CIGAR strings. If you need to guarantee that any random obscure downstream program will be able to handle it. Uptake of CRAM has been pretty slow. Java programs using htsjdk (e.g., picard, IGV and GATK) have only ...


16

The maternal and paternal copies of a chromosome are called haplotypes. Many metazoans (animals) are diploid and have maternal and paternal chromosome contribution during sexual reproduction, not just humans as your question states. Response to Q1 Your question, in other words, is: Why do .bam files not differentiate between haplotypes? Your question ...


15

samtools has a subsampling option: -s FLOAT: Integer part is used to seed the random number generator [0]. Part after the decimal point sets the fraction of templates/pairs to subsample [no subsampling] samtools view -bs 42.1 in.bam > subsampled.bam will subsample 10 percent mapped reads with 42 as the seed for the random number generator.


15

The obvious answer is that different people wrote them. It's fairly common in bioinformatics for people with a computer science background to get frustrated with existing tools and create their own alternative tool (rather than improving an existing tool). Over time, tools with similar initial aims will have popular functionality implemented in them (and ...


15

No, samtools (and therefore bcftools) does not use soft-clipped bases. You can quickly confirm this by using either samtools depth or samtools mpileup to look at a region with a soft-clipped alignment. You'll note that the soft-clipped region isn't used in the depth/pileup (both tools use the same underlying code, so it doesn't matter which you use). If you'...


14

samtools merge merged.bam *.bam is efficient enough since the input files are sorted. You can get a bit faster with sambamba and/or biobambam, but they're not typically already installed and IO quickly becomes a bottleneck anyway.


13

To exclude all possible multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use grep on the uncompressed SAM fields: samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique_mapped.bam Explanation follows... I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM ...


13

Arithmetic on a zero-based coordinate system is less complicated than that on a one-based system, so it appears zero-based is often (not exclusively) used for binary data formats, like BAM or bigBed, or text formats like BED, where computers are used to more efficiently calculate lengths of or set operations on intervals. A more complete answer on the ...


12

You can just set the fields you don't need to *: samtools view -h foo.bam | awk 'BEGIN{FS="\t"; OFS="\t"}{if($1~/^@/) {print $0} else {print "*", $2, $3, $4, $5, $6, $7, $8, $9, "*", "*"}}' | samtools view -bo smaller.bam - This will set each read's name to *, but you can still see where its mate maps with the PNEXT and RNEXT fields. The resulting BAM ...


11

This can be done in R very easily from an indexed .bam file. Given single-end file for sample1. library(GenomicAlignments) library(rtracklayer) ## read in BAM file (use readGAlignmentPairs for paired-end files) gr <- readGAlignments('sample1.bam') ## convert to coverages gr.cov <- coverage(gr) ## export as bigWig export.bw(gr.cov,'sample1.bigwig') ...


11

Arbitrary record access in constant time To get a random record in constant time, it is sufficient to get an arbitrary record in constant time. I have two solutions here: One with tabix and one with grabix. I think the grabix solution is more elegant, but I am keeping the tabix solution below because tabix is a more mature tool than grabix. Thanks to ...


11

Having done 23andme myself I can tell you that your variant file, which contains SNP genotypes, cannot be converted to a bam file. It does not contain the same information as a bam file. It may be helpful to familiarize yourself with those filetypes and the technologies used to obtain them. A SAM/BAM file contains alignments of reads obtained by sequencing. ...


10

Oh you silly sausage, pysam.index takes a bam file name, not a python object. import pysam pysam.index("regular_bwamem_mapping.bam") will index your .bam file.


9

You almost had the correct python code already, you just need to filter out secondary alignments: def get_reads_hist(bam): bam = pysam.AlignmentFile(bam, 'rb') counts = Counter() for query in bam.fetch(): if query.is_secondary: continue nh_count = Counter(dict(query.get_tags())['NH']) counts += nh_count ...


9

The order of -a and -b switched at some point. You want: bedtools coverage -a All_peaks.bed -b file.bam > file.cov.txt For reference, this is the end of the help output in version 2.25: Default Output: After each entry in A, reports: 1) The number of features in B that overlapped the A interval. 2) The number of bases in A that ...


9

use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.


9

samtools quickcheck is all you need. From the manual: Quickly check that input files appear to be intact. Checks that beginning of the file contains a valid header (all formats) containing at least one target sequence and then seeks to the end of the file and checks that an end-of-file (EOF) is present and intact (BAM only). Data in the middle of the file ...


9

It is still slow but grep has a -f option to take in a file samtools view inbam.bam | grep -f read_names.txt > read_locs.txt


9

one-liner Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ($4 >=1) and the 500 in ($4 <=500) to set your window. Change "hg19" to your target sequence. Note, this one-liner does not double-count reads because of uniq. samtools view file_sorted.bam | \ ...


8

Merging sorted files is a linear operation, so any well-implemented tools that do it will do it with approximately the same efficiency. So samtools merge (use the most up-to-date version, as there have been improvements in merge header handling in the 1.3.x and 1.4.x versions), picard MergeSamFiles, etc. These tools need to hold all the input BAM files ...


8

Why not work only with sam files and we need to convert them to bam? SAM is for human only, a binary format like BAM is smaller and parsed much faster by a program (and we want speed and optimize spaces for those huge files, don't we ?). Furthermore, bam files are compressed using a method that allow them to be accessed using random-access, that is how ...


8

The main difficulty here is the use of GRCh38. Unfortunately, despite the fact that it's more than four years old, the major ethnicity-labeled public datasets (1000 Genomes, gnomAD when allele frequencies are enough) still aren't available for that reference. It is necessary to perform a liftover operation, or just use overlapping rsIDs and hope for the ...


8

If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically: samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam | cut -f 5 | tr '[...


8

The quick way to get the number of alignments on each reference is samtools idxstats my_bam.bam Number of reads on each reference is column 3. Although, as has been pointed out, this will give you the total number of alignments per reference, not the total number of reads (each read might give rise to more than one alignment). That said I do tend to us ...


7

It's worth bearing in mind that when outputting compressed BAM, as most tools do by default, they may well be using different levels of compression and/or different libraries, or versions of said, libraries for doing (de)compression which will result in different file sizes. Additionally coordinate sorted BAM will compress more than unsorted BAM. The ...


7

As wkretzsch suggested this was worthy of an actual answer, I feel the obvious solution is missing here; index the FASTQ. Index it As much as I typically hesitate to jump to a solution that requires a script or framework (as opposed to just unix command line tools), there is sadly no samtools fqidx (perhaps there should be), and existing answers suggest a ...


7

GATK has a solution that might work for you: FastaAlternateReferenceMaker, which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records." Input The reference, requested intervals, and any number of variant ROD files. Output A FASTA file representing ...


7

Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your samtools view -H code. With samtools idxstats file.bam you'll get the reads per chromosome (or per nucleotide sequence in your case). For use in featureCounts my advice is to get ...


Only top voted, non community-wiki answers of a minimum length are eligible