Finding allele-specific expressions with AllelicImbalance

An allele-specific expression is a situation that occurs when there is a differential abundance of different allelic variants of a transcript. RNAseq can help to provide quantitative estimates of allele-specific expression for genes with transcribed polymorphisms—that is, variants in the transcript that may result in different proteins. In this recipe, we'll take a look at a method for determining which of the variants of a transcript may have preferential expressions in different samples. The reads will come from different .bam files and the variants must already be known. This implies that you have already carried out a read alignment and a variant call step and have per sample .bam and .vcf files. We'll use the AllelicImbalance and VariantAnnotation packages for this recipe.

Getting ready

You'll need AllelicImbalance and VariantAnnotation from Bioconductor. The AllelicImbalance package provides a small but informative dataset of three SNPs on Chromosome 17 of the hg19 build of the human genome. The files have been extracted into this book's data repository in datasets/ch1/allele_expression .

How to do it...

1. Load libraries and set up an import folder:

library(AllelicImbalance)
library(VariantAnnotation)

region_of_interest <- GRanges(seqnames = c("17"), ranges = IRanges(79478301, 79478361)) 
bam_folder <- file.path(getwd(), "datasets", "ch1", "allele_expression") 

2. Load reads and variants in regions of interest:

reads <- impBamGAL(bam_folder, region_of_interest, verbose = FALSE)

vcf_file <-file.path( getwd(), "datasets", "ch1", "allele_expression","ERP000101.vcf" )
variant_positions <- granges(VariantAnnotation::readVcf(vcf_file), "hg19" )

allele_counts <- getAlleleCounts(reads, variant_positions, verbose=FALSE)

3. Build the ASEset object:

ase.vcf <- ASEsetFromCountList(rowRanges = variant_positions, allele_counts)

reference_sequence <- file.path(getwd(), "datasets", "ch1", "allele_expression", "hg19.chr17.subset.fa")

ref(ase.vcf) <- refAllele(ase.vcf,fasta=reference_sequence)
alt(ase.vcf) <- inferAltAllele(ase.vcf)

3. Run the test on all variants:

binom.test(ase.vcf, n="*")

How it works...

In step 1, the script begins by creating the familar GRanges object describing our region of interest and the folder holding the .bam files of reads.

In step 2, the impBamGAL() function loads in reads in the region of interest. The variant information is loaded into variant_positions—another GRanges object and the reads and variants are used to make allele counts with getAlleleCounts().

With this done, in step 3, we can build the ASESet object, ase.vcf (a class that inherits from RangedSummarizedExperiment), using the constructor function, ASEsetFromCountList(); we then use the setter functions, ref() and alt(), to apply the reference and alternative base identities.

Finally, in step 4, we can apply tests. binom.test() carries out binomial per position per sample (.bam file) tests for deviations from equality in counts in reference and alternative alleles. The parameter n tells the test which strand to consider—in this example, we haven't set up per-strand information, so we use "*" to ignore strandedness.

This will give the following output:

##               chr17_79478287 chr17_79478331 chr17_79478334
## ERR009113.bam          0.500   1.000000e+00   1.000000e+00
## ERR009115.bam         0.125   6.103516e-05   3.051758e-05

There's more...

The preceding analysis can be extended to carry out per strand and per phenotype tests if required. The script would need amending to introduce strand information in the ASESet object construction step. Doing so usually requires that the RNAseq experiment and alignment steps were performed with strandedness in mind and the bioinformatics pipeline up to here configured accordingly. Phenotype information can be added in the construction step using the colData parameter and a vector of phenotype or sample types for columns in the ASESet object.