I’m now on Twitter! @ChimentiMichael
If you like my blog, be sure to follow me on Twitter for updates to this blog, tweets about interesting articles, and brief thoughts about bioinformatics. @ChimentiMichael
Beyond Benjamini-Hochberg: Independent Hypothesis Weighting (IHW) for multiple test correction
Multiple hypothesis testing is a critical part of modern bioinformatic analysis. When testing for significant changes between conditions on many thousands of genes, for instance in an RNA-Seq experiment, the goal is maximize the number of discoveries while controlling the false discoveries.
Typically, this is done by using the Benjamini-Hochberg (BH) procedure, which aims to adjust p-values so that no more than a set fraction (usually 5%) of discoveries are false positives (FDR = 0.05). The BH method is better powered and less stringent than the more strict family-wise error rate (FWER) control, and therefore more appropriate to modern genomics experiments that make thousands of simultaneous comparisons. However, the BH method is still limited by the fact that it uses only p-values to control the FDR, while treating each test as equally powered.
A new method, Independent Hypothesis Weighting (IHW), aims to take advantage of the fact that individual tests may differ in their statistical properties, such as sample size, true effect size, signal-to-noise ratio, or prior probability of being false. For example, in an RNA-Seq experiment, highly-expressed genes may have better signal-to-noise than low-expressed genes.
The IHW method applies weights (a non-negative number between zero and one) to each test in a data-driven way. The input to the method is a vector of p-values (just like BH/FDR) and a vector of continuous or categorical covariates (i.e., any data about each test that is assumed to be independent of the test p-value under the null hypothesis).
From the paper linked above, Table 1 lists possible covariates:
| Application | Covariate |
|---|---|
|
|
|
| Differential expression analysis | Sum of read counts per gene across all samples [12] |
| Genome-wide association study (GWAS) | Minor allele frequency |
| Expression-QTL analysis | Distance between the genetic variant and genomic location of the phenotype |
| ChIP-QTL analysis | Comembership in a topologically associated domain [16] |
| t-test | Overall variance [9] |
| Two-sided tests | Sign of the effect |
| Various applications | Signal quality, sample size |
In simplified form, the IHW method takes the tests and groups them based on the supplied covariate. It then calculates the number of discoveries (rejections of the null hypothesis) using a set of weights. The weights are iterated until the method converges on the optimal weights for each covariate-based group that maximize the overall discoveries. Additional procedures are employed to prevent over-fitting of the data and to make the procedure scale easily to millions of comparisons.
The authors of the method claim that IHW is better powered than BH for making empirical discoveries when working with genomic data. It can be accessed from within Bioconductor.
Unix one-liner to convert VCF to Oncotator format
Here is a handy unix one-liner to process mutect2 output VCF files into the 5 column, tab-separated format required by Oncotator for input (Oncotator is a web-based application that annotates human genomic point mutations and indels with transcripts and consequences). The output of Oncotator is a MAF-formatted file that is compatible with MutSigCV.
#!/bin/bash
FILES='*.vcf.gz'
for file in $FILES
do
zcat $file | grep -v "GL000*" | grep -v "FILTER" | grep "PASS" | cut -d$'\t' -f 1-5 | awk '$3=$2' | awk '$1="chr"$1' > $file.tsv
done
Breaking this down we have:
“zcat $file” : read to stdout each line of a gzipped file
“grep -v “GL000*” : exclude any variant that doesn’t map to a named chromosome
“grep -v “FILTER” : exclude filter header lines
“grep “PASS””: include all lines that pass mutect2 filters
“cut -d$’\t’ -f 1-5” : cut on tabs and keep fields one through five
“awk ‘$3=$2’ : set column 3 equal to column 2, i.e., start and end position are equal
“awk $1=’chr’$1″” : set column one equal to ‘chr’ plus column one (make 1 = chr1)
New paper out: metagenomics study of poultry production environments
I am happy to say that myself and my collaborators in the Department of Occupational and Environmental Health here at the University of Iowa have had our recent work on the bacterial composition of poultry bioaerosols (i.e., the dust that poultry workers breath during their tasks) published in Microbial Biotechnology.
The key figure from this work is the following heat map that illustrates the top taxa that are common to all 21 samples:
What is remarkable about whole-genome shotgun metagenomics is that we are not only surveying bacterial DNA, but also viral, fungal, archaeal, and eukaryotic DNA in one experiment. You can see from the figure that certain viruses are found in all samples, but it is bacteria, particularly Lactobacillus and Salinicoccus, that are the most abundant.
Stay tuned because we will have a paper coming out soon on the fungal composition of these samples as well. In the case of this paper, and our next manuscript, it is the first time whole-genome shotgun metagenomics has been applied to the field of environmental health in poultry environments.