text:
Displaying soft-clipped nucleotides in samtools tview
There are nicer genomics visualization tools available, but the samtools tview command is almost always my go-to for a quick first look at read alignments. I just brought up the following locus in tview.
81055951  81055961  81055971  81055981                                                81055991  81056001  81056011  81056021
GTGGAGGTCGGCAGCAGAGCCTGTGGCAGCT**********************************************CTGAGGGTCGGCGTGGCCTCCTGGTGGAAGTTGCACCTGGTTTGCTT
...............................                                              ...............................................
...............................**********************************************..       ......................................
   ,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,       .....T................................
.....                  ,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
.......                      ,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,                                                                    ...............................................
.............                                                                .........................................A.....
........................                                                     ...............................................
...............................**********************************************..         ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,           ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************.A......                 ,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************..                         ....................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,                         ....................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,    ,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,        ,,,,,
...............................**********************************************..                               ,,,,,,,,,,,,,,
...............................**********************************************..
...............................**********************************************..
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
..................T............**********************************************...............................................
...............................**********************************************...............................................
...............................**********************************************...............................................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************...............................................
...............................**********************************************...............................................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
..............G................**********************************************...............................................
,,c,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
...............................CTAAGCAAGTAACTGTGGTGTTAGAACCTACCTCCCTCGCCGGGCA...............................................
...............................**********************************************...............................................
...............................**********************************************........................................A......
...............................**********************************************...........................................T...
,,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,**********************************************,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
                                                                             ...............................................
                                                                             ...............................................
                                                                             ...............................................
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

At first glance there appears to be only a single read spanning this insertion, although the abrupt end of many reads at the insertion breakpoints made me suspicious. Indeed, upon closer inspection, 19 additional reads have primary alignments with soft clipping at the insertion breakpoints.
By default, soft-clipped portions of the reads are not displayed. Is there any way to toggle this behavior?

How would they be displayed? There's nothing in the CIGAR strings to align the soft-clipped bits to the reference

They could be displayed in the same way that igv does it: as a string of unassigned bases at the end of each read.

A program called gridss does this. It can be found here. It basically extracts the clipped bases and repeatedly realigns them to reference with the aligner of your choice with bwa as default. It will do this until there are no further realignments. The simple command-line for this tool after install is:
java -Xmx512M -cp gridss-VERSION-with-dependencies.jar gridss.SoftClipsToSplitReads I=your_input.bam O=your_output.bam REFERENCE_SEQUENCE=your_reference.fa

I hope this helps!

Could you please explain why you downvoted?

I haven't voted up yet since I haven't had a chance to try out gridss. But yeah, that's an unfair downvote.

This answer does not appear to answer the OP's question?

@winni2k would it be more beneficial for me to have posted an answer that says, "No"? OP's question boiled down to, "Is there any way to toggle this behavior?" There is no way to toggle this behavior, so I proposed an alternate program as a solution. Seems reasonable to me.

Just explaining why someone might down vote. I think it's great that you tried to help though.


text:
Access base aligned to particular reference position
The short version: If I have a SAM record, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?
The long version: I'm using pysam to write some genotyping code. I have a BAM file with reads aligned to some amplicons. I'd now like to genotype a couple of positions of interest in each read while maintaining the phase information. For example, at amplicon 1, if read 1 has A for SNP1 and T at SNP2, then I want to report an A,T genotype for read 1 at amplicon 1.
The typical way to genotype SNPs of interest is to compute a pileup and look at all of the aligned bases position-by-position. The pysam API includes the following example of how to work with pileups.
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    print ("\ncoverage at base %s = %s" %
           (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            print ('\tbase in read %s = %s' %
                  (pileupread.alignment.query_name,
                   pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

However, I lose the phasing information if I analyze the data position-by-position. Is there any way to flip this example so that I'm iterating over reads in the outer loop, and then over positions in the inner loop? In other words, if I'm iterating over reads, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?
As a potential solution, I've though about storing genotype information in a nested dictionary indexed first by read ID and then by position, but depending on coverage that could be incredibly memory intensive.

You can use get_aligned_pairs(); you'll get a list of tuples of (position_in_query_sequence, position_in_reference_sequence)

That's promising! It's been a hot second since I've worked on this code, but it looks like it could indeed allow me to keep phase information on a read-by-read basis. I'll try it out when I can.

You can use the get_reference_positions() function in pysam to get a vector of the positions. Searching that for your positions will allow you to output phased calls.
for read in bam.fetch(some_position):
    positions = read.get_reference_positions():
    if some_other_position in positions:
        idx = positions.index(some_position)
        idx2 = positions.index(some_other_positions)
        ... and so on...

That's not exactly correct syntax, but you get the idea.

Thanks, this is helpful. When there are no insertions (I blocks) in the CIGAR, I can use idx and idx2 to index directly in to record.query_sequence and get the base at the correct position. However, reads with insertions need additional operations to line up the get_reference_positions() indices with the query_sequence indices. So unless I'm missing something, it looks like there's a way to do this without computing a pileup...just not a simple way. :-)

I think you summarized it nicely, possible but not simple.

I had the same problem as you, and so I wrote my own code and sharing it here. No guarantee that it's bug-free ;-)
This code takes advantage of the cigar string. It returns None if the base has been deleted or clipped, or if the bam file did not store the sequence.
from typing import Optional
import pysam

def find_base_in_alignment(alignment: pysam.AlignedSegment,
                           pos: int,
                           bam_stores_revcomp: bool = False) \
        -> Optional[str]:
    idx_q = 0
    idx_r = pos - alignment.reference_start
    if bam_stores_revcomp:
        seq = alignment.query_sequence
    else:
        seq = alignment.get_forward_sequence()

    if seq is None:
        return None

    for op, l in alignment.cigartuples:
        ref_consumed = op in {0, 2, 3, 7, 8}
        query_consumed = op in {0, 1, 4, 7, 8}

        if ref_consumed:
            idx_r -= l
        if query_consumed:
            idx_q += l

        if idx_r < 0:
            if query_consumed:
                # base is in query between idx_q-l , idx_q
                base = seq[idx_q + idx_r - 1]
                return base
            else:
                # position has been deleted
                return None

Let's say you want to genotype all reads by looking at the base at chromosome chrom and position pos:
bamfile = pysam.AlignmentFile(filename, "r" )
for alignment in bamfile.fetch(chrom, pos-1, pos+1):
    base = find_base_in_alignment(alignment, pos)
    # ... do something with it

Note that the parameter bam_stores_revcomp is important. Some aligners will store the original query sequence in the bam, and other aligners (like minimap2) will store the forward sequence (i.e. the reverse complement whenever a read was aligned to the reverse strand). If you use the former, set bam_stores_revcomp to False, otherwise set it to True

Hi, thank you for your code! I think it should be if idx_r <= 0: , then your code would report the exact same data as IGV, excluding insertions. Thank you!

this should be the accepted answer.

This is not [yet] a complete answer, but hopefully it helps get you on the right path.
I wrote some code to retrieve the subsequence from reads in a BAM file that are actually mapped to the reference, which required calculating mapped base locations. Here's a subsection of the code that does this (slightly paraphrased):
  my $refPos = $F[3];
  my $cigar = $F[5];
  $cigar =~ s/[0-9]S$//;
  my $seq = $F[9];
  my $qual = $F[10];
  my $startTrim = 0;
  my $matchLen = 0;
  while($cigar =~ s/^([0-9]+)([MIDNSHP=X])//){
    my $subLen = $1;
    my $op = $2;
    if($op eq "S"){
        $seq = substr($seq, $subLen);
        $qual = substr($qual, $subLen);
      }
    }
    if($op =~ /[M=XI]/){
      $matchLen += $subLen;
    }
  }
  $seq = substr($seq, 0, $matchLen + $addSeq);
  $qual = substr($qual, 0, $matchLen + $addSeq);

So, some things to take note of:

The entire CIGAR string (excluding terminal soft clips) should be parsed up to the required base location
Any soft clips at the start of the sequence need to be accounted for
Only update the read location counter if there's a match, mismatch, or insertion

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Annotating conserved protein domains in bacterial assemblies
I understand the basic principles under which RPS-BLAST operates: search a protein (rpsblast) or nucleotide (rpstblastn) query against a database of protein profiles, a "reverse PSI-BLAST" search. My question has to do with the practical application of this tool to annotate conserved protein domains in bacterial assemblies.
Let's assume I have downloaded NCBI's Conserved Domain Database (CDD) and would like to annotate conserved protein domains in a bacterial genome I have just assembled from Illumina reads. I see two potential approaches.

Run a gene prediction tool on the assembly, translate the predicted gene models, use rpsblast to search the translated sequences against the CDD profiles, and then transform the matches from local protein coordinates to global genome coordinates.
Use rpstblastn to search the genome sequences against the profiles directly.

The latter option would be more compute intensive, since it requires computing a 6-frame translation of each genome sequence. But perhaps it is more sensitive in the presence of small misassemblies?
In practice, is either of these approaches preferred over the other? Are there any considerations I'm missing?

It feels like the arguments would be the same as doing blastx vs. gene prediction + blastp, and i agree with your assessment. perhaps another consideration is do you want to use this rpsblast (using PSSMs/PWMs) or something more complete like profile HMMs

The first option. Using a pre-existing tool, such as PROKKA etc.
Diversity. Bacterial diversity is staggering and species such as E. coli and Salmonella, which are closely related diverged when mammals and dinosaurs/birds diverged, so searching nucleotide sequences is a bad idea as you will get nothing. Especially as you are hoping to find distant sequences that may be missed by the commonly used annotation tools.
Operons. Bacteria, especially thermophiles, often group their protein functionally in operons. Therefore, you can use the principle of guilt by association to guess what a mysterious protein does. It is extremely powerful. Say you have Voet & Voet handy, a pathway with a hole and a mysterious protein you can guess easily what it does.
Note, if you have a protein that you are particularly interested in, I would suggest making your own PSI PSSM file. The automated system cannot read the literature or do some quick checks to see when to stop, so may be overly conservative with the clusters.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Mining documents for nucleotide and peptide sequences
The text mining literature has an emphasis on identifying and normalizing gene names, mutations, pathways, concepts, and so on. I haven't been able to find much, however, on methods for extracting nucleotide and peptide sequences from documents. I did find this tidbit from Aerts et al. (2008), emphasis mine:

Text was split into words and words greater than 10 characters in length with greater than 40% of characters from the capitalized DNA alphabet [ACGT] were extracted using regular expressions to isolate putative DNA sequences.
All putative DNA sequences extracted from each paper were concatenated in the order they appeared in the text into a single fasta sequence and labeled with the corresponding PMID.
Concatenation of sequences was performed to merge sequences split by line breaks in the text conversion, and because we reasoned that inappropriate joins would be reconciled at the genome level by local alignment procedures.
Extracted, concatenated sequences were used as queries to BLAST RepeatMasked versions of genome sequences downloaded from the UCSC genome database...

I'm struck by how simple the author's method for extracting DNA sequences is. But I guess with the limited complexity of the DNA alphabet, that makes sense. This approach wouldn't work with the much more complex protein alphabet.
Are simple regular expressions the text mining community's state of the art for finding DNA/RNA sequences in documents? What methods are used to identify protein sequences?

I am not sure there is a state of the art for that.
Some possible strategies I can think of:

Applying a dictionary to identify language words and discard these.
Using some context information like identifying some words in the same sentence than the sequence.
I guess the format in papers for sequences will not be the same as for text. Maybe you can find a pattern.

I would use an incremental strategy. Use the simplest algorithm then see what it recovers. If there is garbage, find a strategy to remove some garbage with a common pattern. Repeat until you are happy with the result.
Disclaimer: I have some experience in Natural Language Processing.


text:
What does "motif 2 bp or longer" mean?
I'm reading a preprint on a new somatic variant discovery tool, and the first sentence of the Discussion section has me confused.

Across the four datasets analyzed in this study, we discovered that the major source of disagreement between callers originates from somatic variants called within STRs, in particular if the motif is two base pairs or longer.

The "two base pairs or longer" bit is what confuses me. I interpreted "motif" as referring to the short tandem repeats (STRs), which are by definition 2 bp or longer. A motif of < 2 bp is...a homopolymer, not really a motif in the common sense.
Am I missing something here? Or is this simply an awkward word choice on the part of the authors?

I doubt you’re missing anything, I suspect they’re using confused wording.

What's your reference for the definition of STRs? I think it is still ambiguous among the community. Wikipedia states motifs of >=2 bp. However other references include homopolymers as well:
Source 1, source 2 and source 3.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Generic HMM solvers in bioinformatics?
Hidden Markov models (HMMs) are used extensively in bioinformatics, and have been adapted for gene prediction, protein family classification, and a variety of other problems. Indeed, the treatise by Durbin, Eddy and colleagues is one of the defining volumes in this field.
Although the details of each of these different applications of HMMs differ, the core mathematical model remains unchanged, and there are efficient algorithms for computing the probability of the observed sequence given the model, or (perhaps more useful) the most likely hidden sequence given the sequence of observed states.
Accordingly, it seems plausible that there could be a generic software library for solving HMMs. As far as I can tell that's not the case, and most bioinformaticians end up writing HMMs from scratch. Perhaps there's a good reason for this? (Aside from the obvious fact that it's already difficult, nigh impossible, to get funding to build and provide long-term support for open source science software. Academic pressures incentivize building a new tool that you can publish a paper on much more than building on and extending existing tools.)
Do any generic HMM solver libraries exist? If so, would this be tempting enough for bioinformaticians to use rather than writing their own from scratch?

You make it sound like “writing HMMs from scratch” is a huge endeavour. ;-) That said, of course having existing implementations of the relevant algorithms makes sense.

I would also recommend to take a look at pomegranate, a nice Python package for probabilistic graphical models. It includes solvers for HMMs and much more. Under the hood it uses Cythonised code, so it's also quite fast.

There are certainly software libraries for working with HMMs. For a general-purpose implementation in C++, take a look at the SeqAn HMM algorithms.
For your purposes, i.e. “computing … the most likely hidden sequence given the sequence of observed states”, you’d invoke viterbiAlgorithm with your observed sequence and the HMM graph.
More fundamentally I think that most existing, mature implementations are probably found in the domain of signal processing, which has been using them longer than biology, and where most of the underlying theory was developed.

If I remember correctly Ewan Birney's Dynamite (a compiler-compiler) as presented at ISMB 1997 had this functionality, there is also some code here on GitHub https://github.com/birney/wise3 which at least mentions Dynamite.  Suspect Ewan is too busy these days to work on this, although he has tweeted about blowing dust of his old Dynamite, sorry Dynamite code: https://twitter.com/ewanbirney/status/788121636973142016

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Interpreting short indel calls in 1000 Genomes Project VCFs
Consider the following short indel polymorphism rs59679400 on chr7.
$ rsidx search dbSNP/dbSNP_GRCh38.{vcf.gz,rsidx} rs59679400 | cut -f 1-7
[rsidx] running version 0.3
NC_000007.14    12510038        rs59679400      GAAAA   G,GA,GAA,GAAA,GAAAAA    .       .

In the Phase III (2015) 1000 Genomes VCFs (GRCh37), I see the following record for the variant.
$ rsidx search 1000GenomesPhaseIII/ALL.chr7.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.{vcf.gz,rsidx} rs59679400 | cut -f 1-7
[rsidx] running version 0.3
7       12549664        rs59679400      GA      G       100     PASS

That's easy enough to interpret, although it looks like only two of the several possible alleles at that locus are represented.
When I look at the VCFs from the latest NYGC update to the 1000 Genomes data, I see the following. (I'll note that it's really inconvenient that RS numbers are no longer included for any variants.)
$ tabix 1000GenomesNYGC/1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz chr7:12510038-12510038 | cut -f 1-7
chr7    12474780        HGSV_106040     G       <DUP>   .       .
chr7    12510038        7:12510038:GA:G GA      G       .       .
chr7    12510038        7:12510039:GAA:G        GAA     G       .       .

Ignoring the first record, which is a larger structural variant, we see two records at the position in question. My initial question is why they've used two records to encode variation for a single indel variant? As shown in the dbSNP record in the first snippet, VCF allows one to encode an arbitrary number of alternate alleles in a single record.
My second—and more important—question is how one is to interpret these records and resolve the sequence for any particular individual genotyped at the locus? Selecting an arbitrary sample HG00103, I get the following allele calls for the two records? What am I to make of this?
>>> from pysam import VariantFile
>>> vcf = VariantFile("1000GenomesNYGC/1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz")
>>> records = list(vcf.fetch("chr7", 12510037, 12510038))
>>> len(records)
3
>>> records[1].samples["HG00103"].alleles
('G', 'G')
>>> records[2].samples["HG00103"].alleles
('GAA', 'GAA')

this is Marta from the NYGC. Thanks for reaching out with your questions.
The full SNV/INDEL callset that we released (listed as "SNV/INDEL VCFs" in the key resource table in our paper) includes multiallelic representation of all multiallelic sites (i.e. multiple ALT alleles are listed together in column 5 separated by commas) and the record you showed as an example above is listed there in the following way:
chr7    12510038    .   GAAA    GAA,GA,G
The VCF file you were looking at is the phased panel that we generated and released alongside the full unphased callset. The phasing software tools that we used in the paper (SHAPEIT2, SHAPEIT4, Eagle2) as well as most other phasing tools do not support multiallelic sites and require that such sites are split into separate lines and represented as biallelic for the purposes of phasing. As you noted, this is not ideal and can lead to difficulties in interpretation (as others explained above) but it's the best we can do with the existing approaches. We mention this in the Methods section of the paper as well as in the README that's included with the phased panel on the IGSR ftp (all linked in the paper). Please let us know if you have any further questions related to the high-coverage 1kGP resources. Thanks!

Thank you Marta!

Older GATK was using the so-called variant model where overlapping variants were combined during calling. More recent GATK is using the allelic model where one non-reference allele is evaluated at a time.
With the allelic model, a reference allele doesn't mean the sample possesses a reference allele; it means the sample doesn't have the alternate allele. For this particular example, records[1] means a homozygous -A deletion. records[2] means the sample doesn't have the -AA deletion. The two records are consistent with each other.
It is still possible to have conflictive records (for example, homozygous -A deletion and homozygous -AA deletion at the same location). At least this used to be a problem with GATK. Not sure if/how much this has been improved recently.

Often multi-allelic sites would be split into several lines (one per alternate allele) to make filtering easier for specific types of alternate alleles (e.g. filter out MNV vs SNV).
VCF-level, one thing you can do to identify the actual genotype over the locus easier is to join the split alleles into a single multi-allelic record. You can easily do this using bcftools:
bcftools view --samples HG00103 --regions chr7:12510038-12510038 1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz | bcftools norm -m + | bcftools view -H

From here is clear that starting on position chr7:12510038, instead of "GAA", you have homozygous "GA" (1|1), therefore a hmozygous "A" deletion, which could be either the first or second "A" at positions 12510039 and 12510040.


text:
Behavior of `--reference` flag with samtools sort
I'm about to run samtools sort (version 1.6) and I'm reviewing the available configuration options.
$ samtools sort
Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     Set compression level, from 0 (uncompressed) to 9 (best)
  -m INT     Set maximum memory per thread; suffix K/M/G recognized [768M]
  -n         Sort by read name
  -t TAG     Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
  -o FILE    Write final output to FILE rather than standard output
  -T PREFIX  Write temporary files to PREFIX.nnnn.bam
      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
               of OPTION or OPTION=VALUE
  -O, --output-fmt FORMAT[,OPT[=VAL]]...
               Specify output format (SAM, BAM, CRAM)
      --output-fmt-option OPT[=VAL]
               Specify a single output file format option in the form
               of OPTION or OPTION=VALUE
      --reference FILE
               Reference sequence FASTA FILE [null]
  -@, --threads INT
               Number of additional threads to use [0]

I'm curious about the --reference flag. I've never used it before with samtools sort and I'm wondering what behavior it invokes. The online documentation doesn't even mention the flag. Under what circumstances would I want/need to specify the reference genome?

You would only need to supply a reference when using CRAM output encoding. The CRAM format allows storing entries without explicitly encoding the base calls, since one can instead encode the differences to the reference genome.
With BAM output the --reference option has no effect.


text:
Control width of samtools tview "snapshot" when redirecting
I have a large number of loci I would like to examine manually with samtools tview. Rather than typing or copy-n-pasting dozens of coordinates, I was hoping to automate the generation of pileup "images" (text files).
With some experimenting just now, I was happy to learn that if I specify the display option -d H or -d T the output can easily be redirected to a file for later inspection. However, it appears that I cannot control the window captured. Changing the width of the terminal window doesn't change the output, nor does adjusting the -p flag (i.e. -p chr1:501-700). No matter what, tview always reports an 80bp window in text or HTML mode.
Is there any way I can adjust this?

I cross-posted this question on the official samtools issue tracker, and finally got a response.

Looking at the code for the HTML view, I see this and this.
  So maybe try
COLUMNS=200 samtools tview -d H -p CHR:X-Y aln.bam ref.fasta

Partially a shameless plug. See if ASCIIGenome can help you with that.
To review a number of regions that you saved in, say, a bed file, you can either process them in batch or load BAM(s) and bed file and jump from one bed interval to the next using the next -start command.


text:
Index VCF file for rapid access by rsID?
With tabix one can index a VCF file for rapid variant retrieval based on genomic position. I'm wondering if there are any tools that will index a VCF file to allow rapid retrieval using rsIDs and/or other metadata? I'm aware of awk/grep/vcftools one-liners for this purpose, but I'd like to avoid scanning a huge VCF each time I need to retrieve the coordinates of a new batch of rsIDs.

Is re-sorting the VCF files feasible? It's possible to use Tabix on generic tab-delimited files, so sorting and indexing by the RSID column might work.

I attempted this, but it failed on most large VCFs I tried.

John Marshall said:

Moreover note that VCF files are by definition sorted by genomic position, so such an index would just be a hash table of rsID to file offset or genomic position.

I just published rsidx, a new package for indexing VCF files by rsID. As John suggests, I simply created a mapping from rsID to genomic coordinates. These are stored in an sqlite3 database.
It currently takes...several...hours to index build 151 of dbSNP on GRCh38. I'm looking into optimizing sqlite3 parameters to accelerate the indexing, as well as an alternative approach based on a two step process (populate an in-memory sqlite3 db, then dump DB to disk). But in the meantime, the search capability works great.

I am not aware of any widespread index formats for query-by-QNAME for BAM or query-by-ID for VCF. I thought there were historical samtools{-devel,-help} threads to the effect of “in principle one could index a named-sorted BAM file, but there's insufficient demand so no-one's ever implemented it”, but I can't find any just now.
Moreover note that VCF files are by definition sorted by genomic position, so such an index would just be a hash table of rsID to file offset or genomic position.
If you are using dbSNP, probably you have a dbSNP database or API that you can use to query by rsID and get genomic position back. So you can retrieve your VCF contents by rsID in two stages: do that dbSNP query, and use its results to find the variants in your VCF file by position.

high-perf-bio.
I present my high-perf-bio project, which, for the most part, solves the problem of fast extraction of anything values from VCF, BED and unformatted data.
Preparation.
1/ Run the automatic downloading of the toolkit and installation of the dependencies.
wget https://github.com/PlatonB/high-perf-bio/archive/refs/heads/master.zip && unzip -o master.zip && rm master.zip && cd high-perf-bio-master && bash __install_3rd_party.sh

2/ Reboot.
3/ cd path/to/high-perf-bio-master
Create DB with rsIDs.
python create.py -S path/to/gzipped_tables_dir -i rsIDs_field_path

Search rsIDs set in DB.
python annotate.py -S path/to/annotated_gzipped_tables_dir -D db_name -T path/to/trg_dir #-c rsIDs_column_num -f rsIDs_field_path

P.S.1. Detailed help for each high-perf-bio component: python tool_name.py -h
P.S.2. For VCF and BED you can specify the minimum of arguments.
P.S.3. Readme in English is in the early stages of writing. If you have any questions, write to Issues.

This is not too different from Daniel's project - he is more credible on the site plus has written on the  performance of his tool. Adding that information would help us compare your program to his at a glance.

@RamRS
"plus has written on the performance of his tool." The performance of such tools directly depends on the performance of the underlying DBMS. I don't think the speed difference between MongoDB and SQLite is significant for indexed data. "This is not too different from Daniel's project" My toolkit has more functionality. It contains 10 high-level command utilities for solving typical bioinformatics data retrieval tasks. Searching rsIDs is not the only feature of high-perf-bio.

Sure, your tool might be a jack of all trades but for the task at hand, is your tool faster - that is, for the same dataset, does MongoDB give you a significant advantage? In any case, my suggestion is based on the fact that you're not as well knows as Daniel, so providing at least as much info as he does will only help your case.

Take a look at grabix. Using grabix you could implement a binary search on a VCF sorted by rsID and compressed using bgzip.


text:
Load balancing numerous single-threaded processes with Snakemake on the grid
I'm familiar with running Snakemake in a DRMAA-enabled grid environment. (Most often, I write Python programs as wrappers and invoke the workflow via Snakemake's Python API.) Typically, I'm running workflows with rules that perform large, multi-threaded tasks. The rules are usually designed to utilize all available processors on a node, and only a single task is executed on any given node at a time.
But what if I have thousands of single-threaded tasks that I want to distribute among the requested grid nodes? If the corresponding rule(s) are marked threads: 1, is Snakemake smart enough to run multiple tasks on each node? How does it know how many tasks can be run on an individual node—or how does the user specify a limit?


text:
Mix globbing and wildcards when specifying rule input
Consider the following scenario.
├── sample-alice
│   ├── sequence_1.fastq
│   ├── sequence_2.fastq
│   ├── ...
│   └── sequence_n.fastq
└── sample-bob
    ├── sequence_1.fastq
    ├── sequence_2.fastq
    ├── ...
    └── sequence_m.fastq

I'm trying to write a single rule in my Snakefile that will preprocess each sample generically. The rule will execute a single command that takes all of Fastq files associated with the sample simultaneously. The samples do not have the same number of Fastq files.
If it was a single sample, I could do something like this.
rule preprocess:
    input: glob('sample-alice/*.fastq')
    output: 'alice-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

On the other hand, if there was a single Fastq file for each sample, I could use a wildcard to write a generic rule like this.
rule preprocess:
    input: 'sample-{samp}/sequence_1.fastq'
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

Is it possible to combine these two approaches to write a rule for handling multi-file samples generically?

I'm not sure if this exactly solves your question, but in one of my workflows I use the directory as input (as specified in the config.yaml). The shell command then uses {input}/*.fastq.gz See some examples here: https://github.com/wdecoster/nano-snakemake/blob/master/rules/align.smk

A nice solution comes from a colleague of mine, as seen in this Snakemake workflow. The trick is to access the wildcards programmatically using an anonymous (lambda) function. In my example above, it would be implemented as follows.
from glob import glob

rule preprocess:
    input: lambda wildcards: glob('sample-{samp}/*.fastq'.format(samp=wildcards.samp))
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'

Note that your Snakefile needs a from glob import glob for this to work

And it's a bit more succinct using fstrings: glob(f'sample-{wildcards.samp}/*.fastq')


text:
Why does this naive Python paired Fastq reader work?
I wrote some throwaway code for ingesting paired-end reads from a Fastq file into memory.
def ingest(stream1, stream2):
    while True:
        record1, record2 = list(), list()
        for _ in range(4):
            record1.append(next(stream1))
            record2.append(next(stream2))
        yield record1, record2

I had initially put the for loop in a try/Except block, anticipating I'd need to explicitly break out of the endless while loop. However, the code seems to work correctly even after I removed it. I even "corrupted" the final read in one of the paired files, and the ingest function performs reasonably (the last read pair is never yielded).
My question is...why does this code work, since I never explicitly break out of the loop?
Here's an example of how the code is invoked.
>>> with open('seqR1.fastq', 'r') as fh1, open('seqR2.fastq', 'r') as fh2:
...     counter = 0
...     for r1, r2 in interleave(fh1, fh2):
...         counter += 1
>>> counter
5000
>>> r1
['@read5000/1\n', 'AGGTCTGGGGACAGGTGTTTTACAGGCTCCTTTAGGGGTGGCAAGTCTCCCGGAAACATGATTAGTGAAGACACGTTAACAGCGAAACTGACAGAGAACCTCAACATTCTCTGGTAGGAAAAGCCAGGATCCTGGCAGAGGAAGCAACGATTCGAGCAGGGAACTGTGCAGGACAGAGAACAATGCGGCCCCAATTTTGATCCAGGGGTAGTTAGCTTAGCTAAGGAGTTCCCCAGTCCCTAAGTGGAGAGGGGATCCAGAACCCTACAATTTATTATTATATCTTAGCGGATATATCTAT\n', '+\n', 'CCCCCGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG@GGGGGGGFAGGGCGGGGGFGGGGGGGGGCGGGGGDGGGGFGGGGGGGGGCGGGGGGGGGGGGGGGCGGGGGGGGGGGGGGGGG@G<GGGGGGGGGGGGFGGFGFF:DGGFFCGGGGCGGFEGGGGGGEGGGFB:GGGDG=FGGGG,GAGGGDGGFGGGCD3GFG+G?G*FGGG<GCGGGGC@CD<FFGF*:+G=FE<;GG,8C7GF=;9D=G<C5C66G6C>7)FD*6G7D)9;*C:)4D5+9>@G*=0):)*;1*504\n']
>>> r2
['@read5000/2\n', 'TGTGGGGTTCTGGATCCCCTCTCCACTTAGGGACTGGGGAACTCCTTAGCTAAGCTAATTACCCCTGGATCAAAATTGGGGCCGCATTGTTCTCTGTCCTGCACAGTTCCCTGCTCGAATCGTTGCTTCCTCTGCCAGGATCCTGGCTTTTCCTACCAGAGAATGTTGAGGTTCTCTGTCATTTTCGCTGCTAATGTGTCTTCACTAATCGTGTTTCCGGGAGACTTGCCACCCCTAAAGGAGCTTGTTGAACACTTGTCACCAGACCTAAGTCCTTAGGGATGAACCTGTGTCTCGGCAC\n', '+\n', 'CCCCBGC9FGCD-GGGFGGGEFEFF@CCGFGGGGFF+DG6,F69+FGG+FGF<FCC;GGEBFGGEDECFGGG6GDG7BFG9GCG9G9GDGGGG<GDFFFGEGC,<GC=GFG:,:A:EE+GD<:E,FGF=FCAGEFF++G,G,<EGGFGFFA>*FG9EGFCFCGC,B>43GFF@E@2>9ED@,568D=A;C+=5@E/F53D>,6F*9@+,E,C=*0+*@56<*9@:/F,1,:*)*14*0CC=*1)*0*,*)0;+)*))/**)**;0;//9)1+)10.()5/)070*)7*C**3/)(2*****\n']

Your code only runs prior to Python 3.7, now it raises RuntimeError: generator raised StopIteration. This is essentially a bug fix described in PEP 479 -- Change StopIteration handling inside generators. You can enable the new behaviour from Python 3.5 via from __future__ import generator_stop.
Note, it's not related to file-handles per se, you can reproduce it with a smaller example with a simple iterator:
it = iter([1,2])

def f(it):
    while True:
        yield next(it)

print(list(f(it)))

You should not write code relying on the bug

EDIT: Chris_Rands's answer is more authoritative and clarifies that this behavior is essentially a bug.  (I had only tested my example on Python 3.6.8 and hadn't seen PEP 479.)
When next runs out of lines, it raises StopIteration, and apparently that stops the while loop similarly to calling break without you needing to catch the exception.  I couldn't see a definitive source on that but you can see it in action if you catch and re-raise the exception:
def ingest2(stream1, stream2):
    while True:
        try:
            record1, record2 = list(), list()
            for _ in range(4):
                record1.append(next(stream1))
                record2.append(next(stream2))
            yield record1, record2
        except StopIteration:
            print("StopIteration raised")
            raise

You can also try to catch it outside the loop, but it doesn't get re-raised at that point.  (If you try raising one yourself in a while loop, though, it doesn't catch it in the same way, so there must be some extra Python magic going on.  I can't see any difference between the two exception objects so I'd be curious what's going on there.)  Neat example.  You could also make a pair of FASTQ parsers using Biopython's SeqIO.parse for more complex/permanent situations if needed.


text:
How can I output an identity matrix in progressiveMauve?
I'm just getting started with the Mauve aligner and I'm finding the documentation a bit lacking. I'm using the progressiveMauve tool from the command line and would like to output an identity matrix file. The set of output files it generates by default for progressiveMauve --output=outfile are outfile, outfile.backbone and output.bbcols.
The documentation for the tool lists an --input-id-matrix option but this doesn't seem to have any effect. The documentation page for output files does describe an identity matrix file but doesn't say anything about how to generate it.

mauve documentation have a similar command but the usage is --id-matrix=<file>. So very similar to the progressivemauve command (just doesn't include --input). Try it and see if it works.
I think it's a bug for sure as several other people on biostar have complained.

It indeed seems that the --input-id-matrix option doesn't work as expected, as Amar has pointed out.
I tried the following command:
./progressiveMauve genome1.fasta genome2.fasta genome3.fasta --output=test > test.log

And saw the following bits in test.log (truncated):
[...]
Genome conservation distance matrix:
0       0.0767545       0.187454
0.0767545       0       0.193385
0.187454        0.193385        0

[...]
genome content distance matrix:
0       0.0767545       0.187454
0.0767545       0       0.193385
0.187454        0.193385        0

bp distance matrix:
0.9     0.0125  0.075
0.0125  0.9     0.0375
0.075   0.0375  0.9
[...]

Could any of those three tables contain the information you are looking for?
p.s. as an aside, if you are looking into a way to obtain a whole genome similarity matrix you might want to look into mash or sourmash.


text:
Can I ask DESeq2 which variable alone explains which gene's behavior the best?
Imagine I have a dataset of highly periodic gene expression, taken at densely spaced timepoints. I know that the expression is periodic, but different genes have different offsets.
It is trivial to convert my timestamps into a cosine (cos_time = cos(offset + time*pi*2/PERIOD)), so that this new variable is linear with the genes' expressions -- but only with the ones that have a peak at offset.
What if I could have multiple cos_times (let's call them phases) and then ask DESeq2: "OK, which gene's behavior is most linear with which phase?" Barring comparing p-values of each LRT test, is there a better way?
Disclaimer: I know about cosinor, Lomb-Scargle periodograms, ARS, JTK etc, and I know they are likely more applicable in this particular case. I'm interested in whether it is in general possible to ask such a question from DESeq2, even if the variables in question are not periodic. Phase just seemed like the simplest concept to demonstrate this on.

If you have densely spaced time points the using non-linear regression to directly fit the signal to a cosine or sine curve, which would directly produce the offset, would be much simpler. DESeq2 wouldn't be the appropriate tool then, you'd use base R and the nls() function.

If I were to do this, my first guess would be to add in the phases as a separate factor in the statistical model, something like:
~ Condition + sex + phase1 + phase2 + phase3 + ...

And then run LRTs removing phases from the model one at a time? This will leave me with p-values from individual tests, and comparing p-values to each other is probably going to put the wrath of statisticians upon me...
Or is there a way to run DESeq2 in such a way that it spits out the "best phase" for each row?


text:
Single-sample vs. joint genotyping
I am trying to understand the benefits of joint genotyping and would be grateful if someone could provide an argument (ideally mathematically) that would clearly demonstrate the benefit of joint vs. single-sample genotyping.
This is what I've gathered from other resources (Biostars, GATK forums, etc.)

Joint-genotyping helps control FDR because errors from individually genotyped samples are added up, and amplified when merging call-sets (by Heng Li on https://www.biostars.org/p/10926/)

If someone understands this, can you please clarify what is the difference on the overall FDR rate between the two scenarios (again, with an example ideally)

Greater sensitivity for low-frequency variants - By sharing information across all samples, joint calling makes it possible to “rescue” genotype calls at sites where a carrier has low coverage but other samples within the call set have a confident variant at that location. (from https://software.broadinstitute.org/gatk/documentation/article.php?id=4150)

I don't understand how the presence of a confidently called variant at the same locus in another individual can affect the genotyping of an individual with low coverage. Is there some valid argument that allows one to consider reads from another person as evidence of a particular variant in a third person? What are the assumptions for such an argument? What if that person is from a different population with entirely different allele frequencies for that variant?
Having read several of the papers (or method descriptions) that describe the latest haplotype-aware SNP calling methods (HaplotypeCaller, freebayes, Platypus) the overall framework seems to be:

Establish a prior on the allele frequency distribution at a site of interest using one (or combination) of: non-informative prior, population genetics model-based prior like Wright Fisher, prior based on established variation patterns like dbSNP, ExAC, or gnomAD.

Build a list of plausible haplotypes in a region around the locus of interest using local assembly.

Select haplotype with highest likelihood based on prior and reads data and infer the locus genotype accordingly.

At which point(s) in the above procedure can information between samples be shared or pooled? Should one not trust the AFS from a large-scale resource like gnomAD much more than the distribution obtained from other samples that are nominally party of the same "cohort" but may have little to do with each other because of different ancestry, for example?
I really want to understand the justifications and benefits offered by multi-sample genotyping and would appreciate your insights.

Say you are sequencing to 2X coverage. Suppose at a site, sample S has one reference base and one alternate base. It is hard to tell if this is a sequencing error or a heterozygote. Now suppose you have 1000 other samples, all at 2X read depth. One of them has two ALT bases; 10 of them have one REF and one ALT. It is usually improbable that all these samples have the same sequencing error. Then you can assert sample S has a het. Multi-sample calling helps to increase the sensitivity of not so rare SNPs. Note that what matters here is the assumption of error independency. Ancestry only has a tiny indirect effect.
Multi-sample calling penalizes very rare SNPs, in particular singletons. When you care about variants only, this is for good. Naively combining single-sample calls yields a higher error rate. Multi-sample calling also helps variant filtering at a later stage. For example, for a sample sequenced to 30X coverage, you would not know if a site at 45X depth is caused by a potential CNV/mismapping or by statistical fluctuation. When you see 1000 30X samples at 45X depth, you can easily know you are looking at a CNV/systematic mismapping. Multiple samples enhance most statistical signals.
Older methods pool all BAMs when calling variants. This is necessary because a single low-coverage sample does not have enough data to recover hidden INDELs. However, this strategy is not that easy to massively parallelized; adding a new sample triggers re-calling, which is very expensive as well. As we are mostly doing high-coverage sequencing these days, the old problem with INDEL calling does not matter now. GATK has this new single-sample calling pipeline where you combine per-sample gVCFs at a later stage. Such sample combining strategy is perhaps the only sensible solution when you are dealing with 100k samples.
The so-called haplotype based variant calling is a separate question. This type of approach helps to call INDELs, but is not of much relevance to multi-sample calling. Also, of the three variant callers in your question, only GATK (and Scalpel which you have not mentioned) use assembly at large. Freebayes does not. Platypus does but only to a limited extent and does not work well in practice.
I guess what you really want to talk about is imputation based calling. This approach further improves sensitivity with LD. With enough samples, you can measure the LD between two positions. Suppose at position 1000, you see one REF read and no ALT reads; at position 1500, you see one REF read and two ALT reads. You would not call any SNPs at position 1000 even given multiple samples. However, when you know the two positions are strongly linked and the dominant haplotypes are REF-REF and ALT-ALT, you know the sample under investigation is likely to have a missing ALT allele. LD transfers signals across sites and enhances the power to make correct genotyping calls. Nonetheless, as we are mostly doing high-coverage sequencing nowadays, imputation based methods only have a minor effect and are rarely applied.

Thanks, a few followups (broken out over several comments):
Top paragraph sounds more like variant calling than genotyping, which is somewhat infeasible in context of large scale high-coverage studies or ongoing clinical-focused sequencing. Do you see same logic applying to genotyping i.e. conditional on there being a variant at the locus? What I'm struggling with is understanding, batch effect detection aside, why we would trust information on AFS from a somewhat arbitrary collection of samples (genetically speaking) more than a large sample resource like gnomAD?

I don't understand why " Naively combining single-sample calls yields a higher error rate. " Can you elaborate or provide an example, in the context of genotyping?

wrt. to haplotype or assembly-based calling it's simply a reflection of the fact that best methods seem to reason over a region around a locus not just the locus itself, and me trying to figure out where in such a method information from multiple samples can be used. Specifically in the context of  high coverage whole genome sequencing where multiple batches of samples are periodically arriving for analysis and there isn't necessarily a final sample size .

Top paragraph is about deciding the genotype of sample S - it's genotyping. When you know the site allele frequency (AFS is the wrong wording here) in a larger population, you have a better prior. This prior is less accurate across populations, but better than a wright fisher prior. All these theories are only useful to lowCov. For highCov, genotype likelihood has a much bigger effect than prior and cross-sample information.

For the rest, you should ask separate questions. Pooling multiple topics under this one question is hard for you, for me and for readers to comprehend.

Anyway... On assembly-based calling, I already said it has little to do with multi-sample calling. On combining single-sample calls, I already explained in the top paragraph that jointly considering multiple samples helps genotyping – which means not using the info hurts genotyping. BTW, the body of your question is mixing variant calling and genotyping (e.g. FDR is about variant calling), which is causing confusion to me and to Devon below.

The benefit to additional samples is seen in your point 1. The likelihood of making a variant call is a function of (1) the depth of coverage supporting a given variant (ignoring mapping/base quality considerations) and (2) the likelihood of that variant existing given background knowledge. With low depth and no background knowledge, poorly covered variants will be assumed to be sequencing errors. Adding more samples can just serve then to increase the background knowledge on a position.

Thanks Devon, the question specifically targets genotyping i.e. you already have to think that there is a variant there.


text:
Visualize PLINK pedigree files
Is there a pedigree viewer that works with PLINK pedigree files? I saw an old post by Pierre, but I don't see the script on his blog.

We use R package kinship2 here: https://cran.r-project.org/web/packages/kinship2/index.html

We use kinship2 R package to visualize pedigree files (.fam plink files). For some reasons, some characters from plink ped aren't supported (like "-9" for unknown genre or affected status).
Here is the R script we use:
args<-commandArgs(TRUE)
pedPath<-args[1]
pdfFile<-args[2]

library(kinship2)
pedfile <- read.table(pedPath,sep="\t")
pedfile$V5 <- apply(pedfile,1,function(data){if(data[5] == 1 || data[5] ==2){data[5]}else{3}})
pedfile$affected <- apply(pedfile,1,function(data){if(as.numeric(data[6]) == 1){res <- data[6]}else{if(as.numeric(data[6]) == 2){res <- data[6]}else{res <- NA}}})

familyIDs <- as.character(unique(pedfile$V1))

pdf(file=pdfFile)
for(famID in familyIDs){
    ped <- pedigree(id=as.character(pedfile$V2[pedfile$V1==famID]),dadid=as.character(pedfile$V3[pedfile$V1==famID]),momid=as.character(pedfile$V4[pedfile$V1==famID]),sex=as.numeric(pedfile$V5[pedfile$V1==famID]),affected=as.numeric(pedfile$affected[pedfile$V1==famID]))
    plot.pedigree(ped,cex=0.7,align=TRUE)
    title(famID)
}

And we call it this way:
Rscript <kinshipScript.R> <plink ped file> <pdf output>

Now, if only I could convince people to include affected information in the actual .fam files. :-)

@MarkEbbert affected status is not really necessary but it still looks good on the pdf :)


text:
GATK documentation for required depth to reliably call heterozygous mutation in diploid organism?
I'm looking for official GATK documentation (or a recent manuscript) that defines a general recommendation/requirement for sequencing depth to reliably call a heterozygous point mutation in a diploid organism (WGS). In this case, I'm working with humans. Specifically, suppose we're looking at position X in the genome and we want to classify (True/False) whether there was enough coverage to detect a heterozygous point mutation.
What coverage would that be? I'm curious about both: (1) UnifiedGenotyper and (2) HaplotypeCaller
I know general practice is to have 30x coverage genomewide, but I don't see anything in GATK's documentation (i.e., Best Practices) about required coverage.
Really appreciate your help.
Note: I also submitted a question to the GATK team end of last week.

UPDATE: As an update, Sarah Walker (co-author on the poster) responded to my question on the GATK forum. She clarified with the following statement:

We believe the sites around 30X (and above 150X) are being filtered
  due to low mapping quality (since it is whole genome so there are many
  areas that are hard to map), which explains the low sensitivities in
  these areas. As the error bars show, there are very few variants in
  these low mapping quality areas.

She then posted the following plot that falls much more inline with what we would all expect. This is a plot showing sensitivity for NA12878 compared to NIST truth.

shlee, from the GATK forums, pointed me to a poster from AGBT 2018 where the authors compared the HiSeqX and NovaSeq instruments. Figure 3 shows sensitivity for SNPs and short INDELs across multiple allelic fractions (AF) from spike-ins. I was surprised to find that SNPs where AF = 0.4 only had ~50% sensitivity at ~25x coverage, and only ~75% sensitivity at ~50x coverage.
The poster is targeted at somatic variation, but data from figure 3 should behave like a standard diploid germline variant because you often see a wide range of allelic fractions from germline variants. They used spike-ins for figure 3. It's not clear to me whether they used HaplotypeCaller in figure 3, but I don't see how they could have used Mutect2, like they did in figures 1 & 2.
Based on this information, it appears the standard 30x genome-wide coverage is generally insufficient to confidently call many heterozygous mutations, and that ~90-100x is ideal.
Am I missing something? Are we really missing that many heterozygous mutations in standard WGS/WES?
I'll leave this question unanswered for now to allow others to chime in on interpretation, etc. I'd appreciate additional thoughts.


text:
How do you write a .gz fastq file with Biopython?
How do you write a .gz (or .bgz) fastq file using Biopython?
I'd rather avoid a separate system call.
The typical way to write an ASCII .fastq is done as follows:
for record in SeqIO.parse(fasta, "fasta"):
    SeqIO.write(record, fastq, "fastq")

The record is a SeqRecord object, fastq is the file handle, and "fastq" is the requested file format. The file format may be fastq, fasta, etc., but I do not see an option for .gz.
Here is the SeqIO API.

I'm not sure I'm doing it the best way, but here is an example where I read a compressed gzip fastq file and write the records in block gzip fastq:
from Bio import SeqIO, bgzf
# Used to convert the fastq stream into a file handle
from io import StringIO
from gzip import open as gzopen

records = SeqIO.parse(
    # There is actually simpler (thanks @peterjc)
    # StringIO(gzopen("random_10.fastq.gz").read().decode("utf-8")),
    gzopen("random_10.fastq.gz", "rt"),
    format="fastq")

with bgzf.BgzfWriter("test.fastq.bgz", "wb") as outgz:
    SeqIO.write(sequences=records, handle=outgz, format="fastq")

You shouldn't need the StringIO bit for parsing, just make sure you open the gzip handle in text mode under Python3:
records = SeqIO.parse(gzopen("random_10.fastq.gz", "rt"), "fastq")

@peterjc I didn't know about the "text mode", thanks. That explains the error message I got:ValueError: Is this handle in binary mode not text mode?

Here is how I was able to export to a standard gzipped FASTQ file using Biopython. Basically, instead of using SeqIO.write, I directly called the .format method of the SeqRecord object. The example code below imports a gzipped FASTQ file, removes reads that do not contain a G in positions 7, 8, and 9, and writes the results to a gzipped FASTQ file.
from Bio import SeqIO
import gzip

path_in = "path/to/in.fastq.gz"
path_out = "path/to/out.fastq.gz"
handle_in = gzip.open(path_in, "rt")
handle_out = gzip.open(path_out, "wt")

fq = SeqIO.parse(handle_in, "fastq")
for read in fq:
    # Only export reads that have a G in positions 7, 8,
    # and 9
    if read.seq[6:9] == "GGG":
        handle_out.write(read.format("fastq"))

handle_in.close()
handle_out.close()

Using the code example from the original question that converts a FASTA file to a FASTQ file, it would look something like the following:
fastq = gzip.open("path/to/out.fastq.gz", "wt")
for record in SeqIO.parse(fasta, "fasta"):
    fastq.write(record.format("fastq"))
fastq.close()

Caveat: This only works for Python 3. Python 3 differs from Python 2 in how it imports data from binary files (source). In fact, ideally you should be using Python 3.4 or greater since Biopython has compatibility issues with Python 3.3 (source).

As far as I know you can not write compressed fastq on the fly. You will have to write the entire fastq and then compress it.

Hi Aaron Berlin, thanks for your answer and welcome to Bioinformatics stack exchange. Your answer is a bit short; do you have any source you can provide about not being able to write compressed fastq with Biopython, or an alternative python hook that could be used to work around the inability for it to support writing out a compressed file?


text:
Nextflow: cannot get file pairs from Channel.fromFilePairs
I'm remarkably still struggling with some basic Nextflow concepts. I'm simply trying to collect all pairs of .cram and .crai files from a directory using Channel.fromFilePairs(...), but I get no output when I try to view() the Channel.
Question: Why does my Channel.fromFilePairs(...) command not return file pairs, but I get a list of files when using Channel.fromPath(...)?
As far as I can tell, I'm following the same pattern from the following examples:

Nextflow fromFilePairs
.bam/.bai example (except without the extra code)

Code
nextflow.enable.dsl=2

workflow {

    /*
     * Channel.fromPath returns the expected files and prints them to
     * terminal.
     */
    Channel.fromPath('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram')
        | view()

    /*
     * Channel.fromFilePairs does not print anything to terminal.
     */
    Channel.fromFilePairs('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.{cram,crai}')
        | view()
}

Output (note, no file tuples/pairs)
N E X T F L O W  ~  version 21.04.3
Launching `pairs.nf` [pensive_plateau] - revision: f917a5c41b
/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram

TLDR, use:
nextflow.enable.dsl=2

workflow {

    Channel.fromFilePairs('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram{,.crai}')
        | view()
}

Results:
N E X T F L O W  ~  version 20.10.0
Launching `xscript.nf` [special_caravaggio] - revision: a4e53657a0
[A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram.crai]]
[A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram.crai]]

When using the fromFilePairs factory method, the trick is to ensure that the pattern you provide can return a common prefix in order for the files to be grouped as expected. An easy way to check to see if this is working as expected is to directly call the readPrefix helper function. For example:
test_cram = file('/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram')
test_crai = file('/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai')

pattern = '/path/to/adsp/cram/snd10000/A-CUHS-CU000*.{cram,crai}'

println("My CRAM prefix is: ${Channel.readPrefix(test_cram, pattern)}")
println("My CRAI prefix is: ${Channel.readPrefix(test_crai, pattern)}")

You can see that these results do not produce a common prefix:
N E X T F L O W  ~  version 20.10.0
Launching `script.nf` [lonely_jang] - revision: 428df7e55e
My CRAM prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1
My CRAI prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram

But if you change your pattern definition to:
pattern = '/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram{,.crai}'

We can check that it returns a common prefix:
N E X T F L O W  ~  version 20.10.0
Launching `script.nf` [maniac_volta] - revision: a8acf7b03f
My CRAM prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1
My CRAI prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1

Thanks, @Steve. Very helpful. The Channel.readPrefix(...) is a nice test.


text:
Where can I find summary data for how common certain mutation *types* are?
I'd like to know how common certain mutation types are in public data sets like the 1000 Genomes, ExAC, and ESP6500. Specifically, I'd like to know the distribution of stop-gains, stop-losses, frameshift, and other mutation types. For example, what is the median number of stop-gains observed across individuals in the ExAC population? I'd like to have the full distribution, but median values would work.
To be clear, I am not looking for allele frequencies.
I didn't find this kind of summary data when I looked at the respective websites. Is there a good source for these summary data?

Are you looking for statistics on the most severe variant at each genomic location? It is possible for many different transcript variants to map to the same location, which can influence the results that come out of this.

Does Figure 4A help? https://www.nature.com/nature/journal/v491/n7422/fig_tab/nature11632_F4.html

You might be able to calculate yourself the stats you want if you download the data.

For some types of variant annotations, you can pull all the SNPs from dbSNP as they do in this paper.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3500177/

An important clarifying question: Do you only care about this information in the context of humans?

This dbsnp github repo seems to offer that information.
https://github.com/ncbi/dbsnp/tree/master/Build%20Announcements/151


text:
bigWig to bed for regions above/below threshold
I'd like to threshold a bigWig file on a value and get all regions above or below that value in bed format. Ideally without too many file conversions. Are there any tools to do that?

Welcome to Bioinformatics. What have you found? Where did you look up ? (To avoid looking at the same sites) What have you tried? I think that we have some other questions in this site about bed and bigwig. Aren't they helpful?

cross posted: https://www.biostars.org/p/282410/

You can do that with a bit of python:
#!/usr/bin/env python
import pyBigWig
threshold = 10  # Change me
bw = pyBigWig.open("file.bw")  # Change me
of = open("regions.bed", "w")  # Change me

for chrom, len in bw.chroms().items():
    intervals = bw.intervals(chrom)
    for interval in interals:
        if abs(interval[2]) > threshold:
            of.write("{}\t{}\t{}\n".format(chrom, interval[0], interval[1]))
bw.close()
of.close()

You'll need to install pyBigWig (it's available via pip or conda install) and change the lines with Change me in them.

Convert bigWig to WIG via bigWigToWig:
$ bigWigToWig signal.bw signal.wig

Convert WIG to BED via BEDOPS wig2bed:
$ wig2bed < signal.wig > signal.bed

Filter BED by its score column with awk:
$ awk '($5 < 1)' signal.bed > signal.filtered.bed

The test $5 < 1 filters elements from signal.bed that have a score value less than 1. Change this to the threshold test you would want to apply.


text:
How to release an R package with several bed files?
I'm currently creating an R package, and my scripts require before any analysis is done that 1-2 bed files be loaded. Normally, I would run the scripts with the following:
library(data.table)
session_bed = read.table("/path/name/to/bedfile.bed")

and then refer to session_bed throughout the analysis.
However, in order to release this as an R package, I believe that these files must be downloaded at the time on the package installation.
(1) Normally, how are bed files distributed in R packages? Are these gzipped? Where are they placed?
According to http://r-pkgs.had.co.nz/data.html, all .bed files should be inst/extdata I guess. However, it's not clear to me how to "link" to this data within the R package.
How have other packages dealt with this? (Is there a size limit?)
(2) In this case, should the function loading the bed file be executed about loading the R library, library(package_name)? Is there more appropriate option?

You can store it in data/ subfolder. If one of your functions need this bed file, you can include importing it in your function.
Here are some examples how to get the data from data/ subfolder.
There are probably more ways to do this, this is just one way.

As benn said you can store them in data/ subfolder of your package. There is no size limit (AFAIK) but if you plan to publish your package in some repositories they might have limits or other restrictions. Bioconductor requires to divide the data required for a software package at certain threshold, but for small amount of data it is fine.
You can find them via system.file like: read.csv(system.file("data/mydata.csv", package = "mypackage")).
Of course your bed files will need the apropiate function to load them.
You can load when your library is called by modifying the local variable, or using .First (I am not so sure of the last one), but this pollutes the environment and sometimes if preferably to use data function. You might need to convert your files to .Rdata file.


text:
How are centromeric/telomeric regions categorized in RepeatMasker? "Other" category?
I haven't been able to find the following in the documentation at http://www.repeatmasker.org/
(1) Why are certain regions classified as "Other"? Are these regions impossible to concretely classify due to alignment issues?
(2) I believe that centromeric and telomeric regions would be classified as tandem repeats. Are telomeric/centromeric regions classified as "DNA" within the hg19 version of RepeatMasker?
Are these documented somewhere?

I'll leave question 1 to someone else since I don't know.
Regarding question 2, centromeric and telomeric regions are hard masked in hg19, so there's no sequence to repeat mask.
For hg38, the centromeres are labeled ALR/Alpha with a class of centr. Telomeres are still hard-masked, but they tend to have simple repeats next to them (as annotated by repeatmasker, which makes sense).
Apparently the Other class in hg19 is comprised of SVA elements. In hg38, these are reclassified as Retrotransposon class and SVA family, which didn't exist in hg19.

"Apparently the Other class in hg19 is comprised of SVA elements. In hg38, these are reclassified as Retrotransposon class and SVA family, which didn't exist in hg19."

Interesting. Where did you find this information? Aren't SINE elements SVA elements?

Actually, do you have a citation for this? That's probably the most helpful addition to the answer above. "Apparently the Other class in hg19 is comprised of SVA elements. In hg38, these are reclassified as Retrotransposon class and SVA family, which didn't exist in hg19."

I don’t know that there’s anything to cite, this is an observation from looking at the annotations

Ah, you mean you're just looking at the repName for Other and extrapolating these must be SVA elements? You're right---these are SVA_A, SVA_B, SVA_C, SVA_D, SVA_E, SVA_F.

"Hard masked" means we know the sequences but mask them out for analysis purpose. We are unable to assemble most of hg19 centromeres, so we don't know these sequences. I wouldn't use "hard masked" here. We actually don't know hg38 centromeres, either. These are "manmade" sequences. See my answer below.

On (2), centromeres and telomeres primarily consist of satellite sequences. In human centromeres, the dominant type is alpha satellite (aka alphoid). GGAAT micro-satellite (aka SATIII centromeric micro-satellite) is also very common. There are less common types such as beta satellite. Because there are all kinds of centromeric/telomeric sequences of different lengths and properties, RepeatMasker doesn't put them into one class. Note that RepeatMasker has a "centr" class but it does not include all centromeric satellites (no SATIII, for example). In general, it is not advised to find centromeres with RepeatMasker output.
GRCh38 is an unusual reference genome in that its centromeric sequences are computationally generated modeled after Venter as I remember. They are mostly alpha satellite, but are not real sequences (which has always bugged me). GRC know exactly where they put these sequences. GRCh37 is a more typical reference genome. It lacks probably majority of alpha satellites, but still has many. When you map reads to GRCh37, you occasionally see ~10,000-fold read depth even though the average depth is ~30X. These regions are mostly alpha satellites. Mappers just pile all different copies of satellites onto the ones remained in GRCh37.


text:
How to check whether all BAM read contain defined read groups?
I'm trying to investigate whether there are errors within my BAM. After looking at the BAM header to see whether the read groups exists (using samtools, i.e.
samtools view -H file1.bam

Here, it appears that the header includes @RG tags. However, I suspect that the BAM contains reads that either lack the @RG tag or perhaps these read groups are not well-defined, i.e. it looks like some reads in the BAM may be missing an RG tag.
How do I explicitly check whether all reads have an RG tag, and whether these are well-defined?

According to the man page, running samtools stats --split RG <file1.bam> should produce summary statistics separated by read group. If it doesn't produce a list of ungrouped reads, counts/statistics can be compared to running without the --split RG argument.
Here's some example output from a BAM file with combined read groups:
$ samtools stats --split RG both.bam | grep '^SN'
SN  raw total sequences:    2392
SN  filtered sequences: 0
SN  sequences:  2392
SN  is sorted:  1
SN  1st fragments:  2392
SN  last fragments: 0
SN  reads mapped:   2341
SN  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
SN  reads unmapped: 51
SN  reads properly paired:  0   # proper-pair bit set
SN  reads paired:   0   # paired-end technology bit set
SN  reads duplicated:   0   # PCR or optical duplicate bit set
SN  reads MQ0:  0   # mapped and MQ=0
SN  reads QC failed:    0
SN  non-primary alignments: 0
SN  total length:   2292231 # ignores clipping
SN  bases mapped:   2255361 # ignores clipping
SN  bases mapped (cigar):   1538560 # more accurate
SN  bases trimmed:  0
SN  bases duplicated:   0
SN  mismatches: 254720  # from NM fields
SN  error rate: 1.655574e-01    # mismatches / bases mapped (cigar)
SN  average length: 958
SN  maximum length: 9078
SN  average quality:    20.0
SN  insert size average:    0.0
SN  insert size standard deviation: 0.0
SN  inward oriented pairs:  0
SN  outward oriented pairs: 0
SN  pairs with other orientation:   0
SN  pairs on different chromosomes: 0

The --split command has also produced additional .bamstat files, which can be looked at for the summary statistics for individual read groups. In this case I'm only interested in the "raw total sequences" line:
$ grep 'raw total sequences' *.bamstat
both.bam_minusOnly.bamstat:SN   raw total sequences:    1901
both.bam_plusOnly.bamstat:SN    raw total sequences:    491

There are 1901 sequences from the "minusOnly" group, and 491 sequences from the "plusOnly" group, and I know from the previous statistics that there are 2392 total sequences. 1901 + 491 = 2392, so I know there are no ungrouped reads in this file.

I'm not sure how one extrapolates from this whether each of the reads have a RG tag. Could you extrapolate a bit how I would deduce from these statistics?

I've added additional information to demonstrate how this can be done.

using samjdk and invoking the function getReadGroup()

getReadGroup() returns The SAMReadGroupRecord from the SAMFileHeader for this SAMRecord, or null if 1) this record has no RG tag, or 2) the header doesn't contain the read group with the given ID.or 3) this record has no SAMFileHeader

java -jar dist/samjdk.jar -e 'return record.getReadGroup()==null;'  input.bam

the output will be a sam/bam file with missing/bad RG .

A complete solution to validate a SAM/BAM file would be picard's ValidateSamFile.  This will find many other issues with BAM files.
http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile


text:
Extracting the CIGAR string from a BAM via Python?
Is there a standard method in Python to extract a CIGAR string from the BAM?
There are great libraries which parse the CIGAR, e.g. https://pypi.python.org/pypi/cigar/0.1
>>> c = Cigar('10M20S10M')
>>> c.mask_left(10).cigar
'10S20S10M'
>>> c.mask_left(9).cigar
'9S1M20S10M'
>>> Cigar('10S').mask_left(10).cigar
'10S'
>>> Cigar('10H').mask_left(10).cigar
'10H'
>>> Cigar('10H').mask_left(11).cigar
'10H'
>>> Cigar('10H').mask_left(9).cigar
'10H'

It looks like pysam already parses the cigar string, if I'm not mistaken:
import pysam

bam = 'myfile.bam'
bamfile = pysam.AlignmentFile(bam, 'rb')

for read in bamfile:
    if not read.is_unmapped:
        cigar = read.cigar
        print(cigar)

This outputs lists of tuples of the parse output:
[(5, 61), (0, 30), (5, 198)]
[(4, 11), (0, 30), (4, 248)]
[(4, 11), (0, 30), (4, 248)]
....

How could I simply output the CIGAR from read? I would prefer to create a column in a pandas DataFrame

If you really do just want the cigar string then it's read.cigarstring. However, I'm not sure what you're trying to gain with the cigar package from Brent. Unless you want to get the string with the masking changed then the tuple you get from pysam is the same as what you get from cigar (with the exception of the numeric operations instead of character operations).

Oh, the API on the website had me confused with cigar string. Perhaps I could use the parsed cigar from pysam, but I'm confused: How do I parse the list of tuples in order count the total number of insertions, deletions,  matches (and possible other information) from this?  I apologise if this is documented, but the pysam manual is confusing imho.

Should I read the above, e.g. the first line [(5, 61), (0, 30), (5, 198)] as 61 hard-clippings, 30 matches, and 198 hard clippings?

Yes, that is the correct interpretation of that list of tuples, which is documented in the API

I see; this question is largely my misunderstanding of the pysam website. I'm happy to delete this question, unless it's helpful for the community.

# pip install pyranges
# pip install bamread

# wget http://hgdownload-test.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeRikenCage/releaseLatest/wgEncodeRikenCageGm12878CellPapAlnRep1.bam # (half a gig)
# samtools index wgEncodeRikenCageGm12878CellPapAlnRep1.bam

f = "wgEncodeRikenCageGm12878CellPapAlnRep1.bam"

import pyranges as pr
df = pr.read_bam(f, sparse=False, as_df=True)

df[["Chromosome", "Start", "End", "Strand", "Cigar"]].drop_duplicates("Cigar")
#          Chromosome  Start    End Strand             Cigar
# 0              chr1  10579  10606      -               27M
# 36             chr1  29334  29360      +           8M1I18M
# 37             chr1  29334  29360      +           7M1I19M
# 38             chr1  29334  29360      +           5M1I21M
# 44             chr1  29339  29365      +           4M1I22M
# ...             ...    ...    ...    ...               ...
# 19655298       chrM  14550  14575      +       16M5I3M3D3M
# 19660213       chrM  14747  14775      -       15M5D5M4I3M
# 19663987       chrM  15374  15399      -   2M1I2M2D4M3I15M
# 19676864       chrM  16234  16271      +  15M1I4M12D5M1I1M
# 19677003       chrM  16264  16291      +   14M1D3M2D1M3I6M


text:
Downloading an open-source MAF with tumor and normal BAMs?
I'm writing an R package that manipulates MAF files. Normally, I generate a small example of a file and use this for testing. (e.g. a subsetted BAM or MAF).
In this case, I would need a MAF and a tumor/normal BAM pair. My problem is I don't know where I could download an open-source MAF with said BAMs.
Where could researchers download files like this, and not violate privacy laws (i.e. they must be publicly-available BAMs)?

There's the NCI Genomic Data Commons:

For all submitted sequence data, including BAM alignment files, the GDC generates new alignments in BAM format using the latest human reference genome GRCh38 with standard alignment pipelines. Using these standard alignments, the GDC generates high level derived data, including normal and tumor variant and mutation calls in VCF and MAF formats, and gene and miRNA expression and splice junction quantification data in TSV formats.

But that needs an NIH account for individual-level data, which might be too restrictive for you.
There's also the Texas Cancer Research Biobank (7 human cancer patients with matched normal):

These tumor and normal specimen sequence files (FASTQ and BAM), somatic variant calls (VCF and MAF) and germline MAF files are available for each consented patient as described on the TCRB Open Access Privacy page.

That also requires account registration, but the registration process is automated and much simpler (requires just a name, and acceptance of usage conditions). There's a published paper about that dataset.
I had some vague feeling that another repository existed (created within the last year), but I can't seem to find it from quick web searches.


text:
How to concatenate "by chromosome"-VCFs?
I have a several VCFs which are VCF which only contain information by chromosome. That is, there's a chromosome 1 VCF (with only chr1), a chromosome 2 VCF (with only chr2), etc.
I checked to make sure that these VCFs were valid via VCFtools, i.e.
$ vcf-validator chr1.vcf

which works---these are valid VCFs I was given.
Now, I would like to combine these VCFs into one VCF.
I naïvely tried the following cat operation:
$ cat chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf > total_chroms.vcf

This doesn't work properly though
$ vcf-validator total_chroms.vcf
The header tag 'contig' not present for CHROM=chr1. (Not required but highly recommended.)
Could not parse the line, wrong number of columns: [##fileformat=VCFv4.2\n]
 at /path/vcftools-0.1.14/perl/Vcf.pm line 172, <__ANONIO__> line 191016.
     Vcf::throw('Vcf4_2=HASH(0x1ae7208)', 'Could not parse the line, wrong number of columns: [##filefor...') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 335
     Vcf::next_data_array('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 3457
     Vcf4_1::next_data_array('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 2574
     VcfReader::run_validation('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14//bin/vcf-validator line 60
     main::do_validation('HASH(0x16ada68)') called at  /path/vcftools-0.1.14//bin/vcf-validator line 14
$

What tools are available to merge these VCF's together into a total VCF?

I would recommend bcftools concat.  You can't just cat them together because each file has a header section.  The bcftools command will handle all that for you.  Each vcf file must be sorted prior to calling concat
bcftools concat -o total_chroms.vcf chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf

It's a bit odd. vcf-validator chr1.vcf shows no problems. I believe the output of bcftools concat -o chroms12.vcf chr1.vcf chr2.vcf shows an error: FIXME: sequence name chr1 in chr1.vcf

@ShanZhengYang the errors seem to indicate that you have chr1 in the CHROM field for some rows, but no corresponding line in the header. Can you check if that is true in the chr1.vcf file? Maybe the validator somehow is missing this error?

@juod The VCF header looks correct to me. #CHROM POS ID REF ALT QUAL FILTER INFO

@juod "but no corresponding line in the header" Maybe I don't understand this. Besides the code pasted above, there is no other "chrom" information in the header

According to the VCF spec https://samtools.github.io/hts-specs/VCFv4.1.pdf there should be header lines for each chromosome this is from the first example in the linked documents - ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>

@Bioathlete I see. I think I understand now---this issue has been discussed before: https://github.com/samtools/bcftools/issues/326
I know that these VCFs are from hg38. Could I simply add the line ##contig=<ID=chr1,length=248956422>, etc. to the VCFs in question? I would just take the chromosome lengths from hg38

@ShanZhengYang try that - it should work. I don't even think it's necessary to provide the length, but you could experiment and let us know.

@juod The above worked! Thanks for the help!

use picard GatherVcf: http://broadinstitute.github.io/picard/command-line-overview.html#GatherVcfs

Gathers multiple VCF files from a scatter operation into a single VCF file. Input files must be supplied in genomic order and must not have events at overlapping positions.

The error I get here is Exception in thread "main" picard.PicardException: In order to index the resulting VCF input VCFs must contain ##contig lines.. It is true, my "by chromosome"-VCFs do not contain a ##contig line in the header...but I'm not sure that would matter

Edit above: Apparently this issue has been discussed before: github.com/samtools/bcftools/issues/326
Question: how could I add these contig lines if I don't know this information?I know that these VCFs are from hg38. Could I simply add the line ##contig=<ID=chr1,length=248956422>, etc. to the VCFs in question? I would just take the chromosome lengths from hg38

"but I'm not sure that would matter" it checks the order of the contigs using the dict. use UpdateVcfSequenceDictionary https://broadinstitute.github.io/picard/command-line-overview.html

"it checks the order of the contigs using the dict. use UpdateVcfSequenceDictionary". Ah, I understand now. Thanks for the help!

The best tool for the job is probably bcftools as suggested by Bioathlete, but you can also do this manually. You just need to collect all header lines from all vcf files, remove any duplicates, and then print all the headers + the actual data into the new one:
grep '^##' chr*vcf | sort | uniq > all.vcf
grep -m1 '^#CHR' chr1.vcf >> all.vcf  ## Get the chr header line
grep -v '^#' chr*vcf >> all.vcf

I appear to get a malformed VCF header with the above commands

@ShanZhengYang malformed how? Could you edit your question and show us an example of your files?

So, I see the above grep commands give the same behavior with all of my standard VCFs. I'll get something like this in all.vcf:
chr10.vcf:##contig=<ID=chr10,length=133797422> chr10.vcf:##fileDate=20121011 chr10.vcf:##fileformat=VCFv4.2 chr10.vcf:##reference=path/reference.fa chr10.vcf:##source=foobar chr11.vcf:##contig=<ID=chr10,length=133797422> chr11.vcf:##fileDate=20121011 chr11.vcf:##fileformat=VCFv4.2 chr11.vcf:##reference=path/reference.fa chr11.vcf:##source=foobar ....


text:
R package development: How does one automatically install Bioconductor packages upon package installation?
I have an R package on github which uses multiple Bioconductor dependencies, 'myPackage'
If I include CRAN packages in the DESCRIPTION via Depends:, the packages will automatically install upon installation via devtools, i.e.
devtools::install_github('repoName/myPackage')

This is discussed in Section 1.1.3 Package Dependencies, in Writing R Extensions
Is there a way to streamline this such that packages from Bioconductor are automatically installed as well?
Normally, users install Bioconductor packages via BiocLite, e.g.
source("http://www.bioconductor.org/biocLite.R")
biocLite("edgeR")

Note that as of Bioconductor release 3.8 onwards, bioconductor packages are installed with BiocManager::install() which can be installed from CRAN

There's a trick to this where one needs to add biocViews: to the package Description. That's the only solution I've ever seen to allowing automatic installation of bioconductor dependencies. If you need a couple of examples, then click through the link I posted and scroll down to pull requests referencing that issue, they will generally include the actual example.

I didn't find the examples in the pull requests you mentioned. However, I did search the github mirror on CRAN. There are DESCRIPTION files with a lone biocViews: line before Depends:. Others include tags, like biocViews: Infrastructure, RNA

Maybe we should post an example DESCRIPTION so it's absolutely clear to future biohackers how this is done...this news needs to spread!

You don't need anything after biocViews:, that alone suffices (see ?devtools::package_deps).

As suggested, here’s an example showing the relevant lines from a DESCRIPTION file from a CRAN/GitHub hosted project that has Bioconductor dependencies (truncated):
Depends:
    R (>= 3.3.0)
biocViews:
Imports:
    methods,
    snpStats,
    dplyr

The relevant bit is the empty biocViews: declaration, which allows the Bioconductor dependency {snpStats} to be automatically installed.

r Check() would still fail for me with "Packages required but not available"

@Edmondo1984 Hard to know what you’re doing wrong without seeing your code. As evidenced by the linked package (which can run R CMD check without failure), this works just fine.

I have a clear repro. I have a DESCRIPTION containing biocViews: and an Import on a bioconductor package. I do a check() it fails, I install the package manually and do a check, it works. The problem is, how is this supposed to work on a CI/CD server? Shouldn't R download the packages itself on a check()?

@Edmondo1984 No, R CMD check doesn’t install packages. The same is true for CRAN packages. If you want to run check on CI, the conventional way is to first install the package and its dependencies, and then run check on it.

Thanks, what is the right way to install the package and its dependencies?

@Edmondo1984 The easiest way is probably via remotes::install_deps(dependencies = TRUE). Check out r-lib/actions, it’s a comprehensive repository of CI/CD actions for R on GitHub. Even if you’re not using GitHub Actions, these can serve as a template for you.

I now recommend using pak::pak() instead of remotes::install_deps().

This appears to only work for remotes::install_github(). If I try install.packages('PhenotypeSimulator') there will be a warning: Warning in install.packages: dependency ‘snpStats’ is not available.

@Axeman Yes. But as mentioned in my previous comment, the recommended solution now is pak::pak(). The ‘remotes’ package is for all intents and purposes obsolete and superseded by ‘pak’.

I have a counter example to other answers. BiocViews is not needed to install Bioconductor dependencies on CRAN. I have a package on CRAN that suggests Bioconductor packages and they are found without using BiocViews.
If you are installing a package from github and it has dependencies on Bioconductor you can use directly BiocManager::install("user/package") and it will use remotes::install_github internally and pull also packages from Bioconductor (from the right release).
BiocViews is/was used to detect that a package is for Bioconductor by remotes, devtools and sessioninfo among other packages that help the development and testing of packages.
In addition Bioconductor doesn't add a Repository: Bioconductor on the provided DESCRIPTION field of the packages, while CRAN adds Repository: CRAN.

Alternatively you also can set options(repos=BiocManager::repositories()) (e.g. in your .Rprofile which will make Bioconductor behave like a CRAN repository (there actually are multiple CRAN-style repositories in Bioconductor) after which you can simply run install.packages() for all packages irrespective if it is hosted on CRAN or Bioconductor.
As a consequence, packages with bioconductor dependencies will automatically resolve those dependencies without any of the aforementioned additions.


text:
Is there a standard method to validate whether a bigWig file is correct?
I have a bigWig file of extension *.bw (which is standard), but I am suspicious that this file may be of an incorrect format.
Is there a standard method in the field to validate a *.bigWig file?

As noted in another answer, it's relatively easy to tell if a file is bigWig or not by opening it up with, say, standard UCSC Kent utilities like bigWigToBedGraph or bigWigToWig, or other bigWig readers, and seeing if you get data back.
Depending on what you're doing with the bigWig file, the trickier part can be how the original data was indexed.
If the bigWig file was originally made from a bedGraph or BAM file, its coordinates will be 0-indexed when converted to human-readable form. If the bigWig file was originally made from a wig file, its coordinates will be 1-indexed when converted to a human readable form. How the file's coordinates were originally indexed may change how you work with them downstream for set operations, for instance.
I think the only way to know about the source is to look at the header contained within the bigWig file, after extraction to human-readable form with UCSC tools. A bigWig file made from a bedGraph file will generally contain the header #bedGraph section <coords> in it somewhere (after conversion), for example.

You can actually have multiple representations (bedGraph and wiggle) within a single bigWig file, provided they're in different blocks. If Kent's tools don't enforce 0-based coordinates then I would consider that a bug, as the bigWig structure itself is 0-based.

There's no standard validation method, but a simple one would be python -c "import pyBigWig; bw = pyBigWig.open('your_file.bw'); print(bw.IsBigWig());". If you get True then it's a bigWig file. If you get False or an error then it's not (False indicates a bigBed file).
Alternatively, file your_file.bw should return data.

Link for future readers: https://github.com/deeptools/pyBigWig