knitr::opts_chunk$set(echo = TRUE)
The variant and genotype calling pipeline described in this vignette is intended for recent or ancient allopolyploid species, in which the reference genome sequence contains many paralogous regions that do not recombine with each other at meiosis. In highly duplicated genomes such as these, conventional alignment software, especially when restricted to only returning one alignment, will often align sequences to the incorrect paralog, causing issues with variant and genotype calling downstream. The pipeline described in this vignette uses read depth distribution across a population (the $H_{ind}/H_E$ statistic) to determine whether a group of tags is behaving like a Mendelian locus. Attempting to optimize this statistic, it rearranges tags among paralogous alignment locations to correct alignments and generate well-behaved markers.
This pipeline is currently only designed for natural populations, diversity
panels, or mapping populations in which the most recent generation was produced
by random intermating among all progeny. F1, self-fertilized, and backcrossed
populations will not be processed correctly because expected heterozygosity
cannot be estimated directly from allele frequencies in these populations.
However, the $H_{ind}/H_E$ statistic can still be used from within polyRAD
to filter non-Mendelian markers from these populations using the HindHeMapping
function.
If you use this pipeline, please cite:
Clark LV, Mays W, Lipka AE, and Sacks EJ (2020) A population-level statistic for assessing Mendelian behavior of genotyping-by-sequencing data from highly duplicated genomes. bioRxiv doi: 10.1101/2020.01.11.902890, Version 1.
Below is an overview of the pipeline. It uses a combination of existing bioinformatics tools, custom Python scripts included with polyRAD, and R functions within polyRAD.
process_sam_multi.py
to group tags based on sets of alignment
locations.readProcessSamMulti
. Run HindHe
on the imported data to identify
problematic samples and estimate inbreeding.process_isoloci.py
to sort tags into their correct alignment locations
within the groups of paralogs identified by process_sam_multi.py
, and filter
out tags that cannot be adequately sorted.readProcessIsoloci
, then perform genotype calling.Everything except for the R portion of this pipeline will need to be run on your operating system's Terminal/Shell/Command Prompt. Bowtie2 requires Mac or Linux. Python 3 is needed for the Python portions of the pipeline. If any of this presents a barrier to you, I recommend finding someone experienced with bioinformatics at your institution who can give you some advice.
Line breaks have been inserted in some terminal commands for readability in this tutorial. (Those same line breaks should not be typed at the command line.)
An overview of the TASSEL-GBSv2 pipeline is available at https://bitbucket.org/tasseladmin/tassel-5-source/wiki/Tassel5GBSv2Pipeline. TASSEL can be downloaded from https://www.maizegenetics.net/tassel.
We will only run the first couple steps of the TASSEL-GBS pipeline. The files that you will need to start with are the key file and your FASTQ files.
-c
parameter to something higher than the default
(10 is still pretty conservative but will discard most sequencing errors) to
save some processing downstream.The commands might look like (on Windows):
run_pipeline.bat -fork1 -GBSSeqToTagDBPlugin -e PstI-MspI -i D:\Msa\raw_data\ -db D:\Msa\Msa.db -k D:\Msa\key.txt -kmerLength 80 -endPlugin -runfork1 run_pipeline.bat -fork1 -TagExportToFastqPlugin -c 10 -db D:\Msa\Msa.db -o D:\Msa\Msa_tags.fq -endPlugin -runfork1 run_pipeline.bat -fork1 -GetTagTaxaDistFromDBPlugin -db D:\Msa\Msa.db -o D:\Msa\Msa_ttd.txt -endPlugin -runfork1
You may need to adjust the memory parameters to accomodate a large dataset.
If for some reason you can't use TASSEL-GBS (for example, your FASTQ files don't have inline barcodes, or your protocol doesn't involve restriction enzymes) you might be able to craft a custom alternative for yourself. What you will need out of this step is a file formatted identically to a TagTaxaDist file from TASSEL, and a FASTA or FASTQ file containing the same sequence tags as the TagTaxaDist file.
A TagTaxaDist file is a tab-delimited text file formatted as below, and contains the read depth of every tag in every sample (taxon). I have used a simple convention for sample names here, but the names can be anything you want. Depending on the complexity of your dataset, it can be quite a large file.
| Tag | Sam1 | Sam2 | Sam3| Sam4 | |-----|------|------|-----|------| | TGCAGAAATCATAGATTAAGGATAT| 0 | 24 | 1 | 0 | | TGCAGAACAGGATACGATACCCCTT| 2 | 0 | 0 | 0 | | TGCAGAACAGTATACGATACCCCTT| 0 | 5 | 0 | 0 |
Next, align the FASTQ file generated in step 2 of the previous section to your reference genome. If you have never used Bowtie2 on your reference genome before, you will first need to make an index.
bowtie2-build Msinensis_497_v7.0.hardmasked.fa Msi_DH1_7
Importantly, you will need to set the -k
argument to something higher than the
number of subgenomes within your reference genome. For example, in an
allotetraploid do at least -k 3
, and in an allohexaploid at least -k 4
.
For my Miscanthus dataset, I ran:
bowtie2 -k 3 --very-sensitive -x Msi_DH1_7 -U Msa_tags.fq -S Msa_align.sam
Next you will run process_sam_multi.py
. To locate this file, in R run
system.file("python", "process_sam_multi.py", package = "polyRAD")
For convenience, you might copy this file to somewhere closer to your working
directory where you have the files generated in the previous step. While you
are at it, also copy isoloci_fun.py
and process_isoloci.py
to somewhere
convenient.
As input, you will need the SAM file generated by Bowtie2 and the TagTaxaDist
file from TASSEL. Be sure to set the -g
argument to the number of subgenomes
if it is different from the default of 2 (allotetraploids); for example, in
an allohexaploid you should set -g 3
.
For descriptions of other arguments, see
python process_sam_multi.py --help
The -c
argument is useful if you want to divide your data into smaller chunks
for downstream processing (if you are limited on RAM or have many processors at
your disposal). If there were some samples in your TASSEL database that you
already know you want to get rid of, you can list the samples that you want to
keep in a text file (one sample name per line) and pass that to the -s
argument.
In my Miscanthus data, the command looked like:
python process_sam_multi.py Msa_align.sam Msa_ttd.txt split_depths/Msa_split -c 5
The output will be a pair of CSVs for each chunk. One is named "align", containing the set of alignment locations for each tag, with their respective number of mutations and CIGAR strings. The other, named "depth", is a subsetted version of the TagTaxaDist file, containing the same tags in the same order as "align".
Before proceeding, we will want to get a look at the data to identify any outlier samples that may be a different ploidy from the rest, interspecific hybrids, or highly contaminated. I recommend removing them to avoid biasing $H_{ind}/H_E$ estimates, but keep in mind that this means they will be excluded from all downstream analysis. If a sample is important to your study and you know of a reason why it would have a different heterozygosity from most samples, you should keep it.
We will import just 1000 loci using readProcessSamMulti
. Sequence tags will
be assigned preliminary alignment locations based on where they had highest
sequence similarity, or to a random location if there was a tie. This is
essentially a preview of what the data would look like if we went ahead with
variant calling using the top alignments returned by Bowtie2.
library(polyRAD)
myRADprelim <- readProcessSamMulti("Msa_split_1_align.csv")
We will then estimate a $H_{ind}/H_E$ matrix for this data, as well as getting a sum of the read depth for each individual.
# subset the object to have diploids and a few tetras diploids <- readLines("diploids.txt") myRADprelim <- SubsetByTaxon(myRADprelim, c(diploids, "KMS397", "KMS444", "UI11-00032"))
hh <- HindHe(myRADprelim) TotDepthT <- rowSums(myRADprelim$locDepth)
load(system.file("extdata", "MsaHindHe.RData", package = "polyRAD"))
Now we will look at the distribution of $H_{ind}/H_E$ across samples, and remove outlier samples.
hhByInd <- rowMeans(hh, na.rm = TRUE) plot(TotDepthT, hhByInd, xlog = TRUE, xlab = "Depth", ylab = "Hind/He", main = "Samples") abline(h = 0.5, lty = 2)
For a diploid population, we aren't expecting values above 0.5. (In general, you shouldn't see values above $\frac{ploidy - 1}{ploidy}$.) So, let's check which are the outlier samples, then remove them from the dataset.
threshold <- mean(hhByInd) + 3 * sd(hhByInd) threshold hhByInd[hhByInd > threshold] hh <- hh[hhByInd <= threshold,]
myRADprelim <- SubsetByTaxon(myRADprelim, rownames(hh))
If you remove any samples like this, you should also make a text file indicating samples to retain in the analysis, which we will use in the next Python step.
writeLines(rownames(hh), con = "samples.txt")
Now we can get a sense of what $H_{ind}/H_E$ looks like across loci. From this
we can estimate inbreeding and get a rough sense of what proportion of loci
will need to be adjusted or filtered by process_isoloci.py
.
hhByLoc <- colMeans(hh, na.rm = TRUE) hist(hhByLoc, breaks = 50, xlab = "Hind/He", main = "Loci", col = "lightgrey")
Again, we expect values below 0.5 (i.e. $\frac{ploidy - 1}{ploidy}$). The tail that we observe above that value likely represents groups of tags that aligned to the same location but in fact represent paralogous loci. The values near zero likely represent loci with very high amplification bias or overdispersion. We will take the peak at 0.3 to represent typical well-behaved markers.
To estimate inbreeding, we'll also need to estimate overdispersion. You may want to check a broader range of values than 9 to 14. In the matrix of quantiles, we are looking for columns where the 25th percentile is about 0.25, the 50th percentile is about 0.5, etc.
overdispersionP <- TestOverdispersion(myRADprelim, to_test = 9:14) sapply(overdispersionP[names(overdispersionP) != "optimal"], quantile, probs = c(0.01, 0.25, 0.5, 0.75, 0.99))
cat("Optimal value is 12.\n", sep = "\n") print(structure(c(0.025793146160144, 0.253340971566736, 0.552888349753057, 0.84390807760186, 1, 0.019449092419997, 0.241258741258741, 0.541687845916098, 0.836394612538573, 1, 0.0147776703743338, 0.226466237547336, 0.523304087800422, 0.829407211200752, 1, 0.0114170501660437, 0.211650502292279, 0.505828679120272, 0.82350760103819, 1, 0.00893782083692289, 0.197336617073809, 0.490038344650588, 0.818544664007547, 1, 0.00697809321403532, 0.185307346326836, 0.4777134652488, 0.813849077168056, 1), .Dim = 5:6, .Dimnames = list( c("1%", "25%", "50%", "75%", "99%"), c("9", "10", "11", "12", "13", "14"))))
If $H_{ind}/H_E$ is 0.3 and overdispersion is 12, we can estimate inbreeding to be about 0.4 using the figure below, which was generated using simulated data.
If all markers were Mendelian, we can see what the expected distribution of $H_{ind}/H_E$ would look like, given the sample size and read depth distribution observed in the dataset.
ExpectedHindHe(myRADprelim, inbreeding = 0.4, ploidy = 2, overdispersion = 12)
message("Simulating rep 1") message("Completed 5 simulation reps") cat(c("Mean Hind/He: 0.286", "Standard deviation: 0.0891", "95% of observations are between 0.144 and 0.481"), sep = "\n") load(system.file("extdata", "MsaHindHe3.RData", package = "polyRAD")) hist(testhhdist, xlab = "Hind/He", main = "Expected distribution of Hind/He", breaks = 30)
We are now ready to run the script process_isoloci.py
. It will be run
individually on each pair of "align" and "depth" files generated by
process_sam_multi.py
. Note that isoloci_fun.py
must be present in the same
directory as process_isoloci.py
. We'll use an expected $H_{ind}/H_E$ of
0.3 and a maximum of 0.5, rounded from the output of ExpectedHindHe
above.
python process_isoloci.py split_depths/Msa_split_1_align.csv split_depths/Msa_split_1_depth.csv -e 0.3 -m 0.5 -s samples.txt
This will take some time to process, and will output a file called
Msa_split_1_sorted.csv
in the split_depths
folder. This file will contain
corrected alignment results, tag sequences, position and sequence for the
variable portion of the tag with respect to the reference genome, and depth
data for all retained samples. Tags are arranged into isoloci (i.e.
assigned to alignment locations) to optimize $H_{ind}/H_E$ first and number of
mutations from the genome second. Isoloci are removed from the dataset if
they exceed the expected $H_{ind}/H_E$ by a certain amount (essentially, if
it is closer to what would be expected for double the ploidy than for the
input ploidy).
We can now import the sorted dataset into polyRAD. Since there is a small
number of individuals in this dataset, we will lower the filtering thresholds
for loci we import. The possiblePloidies
argument is left at the default of
2
because although the species is allotetraploid, the sorted loci should now
behave in a diploid fashion.
myRAD <- readProcessIsoloci("Msa_split_1_sorted.csv", min.ind.with.reads = 80, min.ind.with.minor.allele = 5)
We can see what the distribution of $H_{ind}/H_E$ looks like after we have sorted and filtered isoloci.
hh2 <- HindHe(myRAD) hh2ByInd <- rowMeans(hh2, na.rm = TRUE) hh2ByLoc <- colMeans(hh2, na.rm = TRUE)
load(system.file("extdata", "MsaHindHe2.RData", package = "polyRAD"))
hist(hh2ByInd, xlab = "Hind/He", main = "Samples", breaks = 20, col = "lightgrey") hist(hh2ByLoc, xlab = "Hind/He", main = "Loci", breaks = 50, col = "lightgrey")
By default, readProcessIsoloci
runs MergeRareHaplotypes
internally. (This
behavior can be turned off with mergeRareHap = FALSE
.) Those
loci that now appear to have a very inflated $H_{ind}/H_E$ are ones where all
true minor alleles were rare enough to be merged into the common allele, and
any remaining minor alleles only represent sequencing error or alleles that were
amplified very poorly. We can get rid of these.
mean(hh2ByLoc <= 0.5) # proportion of loci retained keeploci <- names(hh2ByLoc)[hh2ByLoc <= 0.5] head(keeploci) hist(hh2ByLoc[keeploci], xlab = "Hind/He", main = "Loci", breaks = 50, col = "lightgrey")
myRAD <- SubsetByLocus(myRAD, keeploci)
Genotype calling can then be performed as normal.
myRAD <- IteratePopStruct(myRAD, overdispersion = 12)
If you want genotypes output as traditional SNP markers, you can use the
RADdata2VCF
function.
RADdata2VCF(myRAD, file = "Msa_test.vcf")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.