Tools for working with amplicon sequencing data.
Developed and maintained by Ben J. G. Sutherland, Ph.D. (Sutherland Bioinformatics & Vancouver Island University).
The development of this pipeline has been supported by the following organizations: Support and Funding page.
Note: this software is provided 'as is', without warranty of any kind, express or implied, including but not limited to the warranties of merchantability, fitness for a particular purpose and noninfringement. In no event shall the authors or copyright holders be liable for any claim, damages, or other liability, whether in action of contract, tort or otherwise, arising from, out of, or in connection with the software or the use or other dealings in the software.
- Ion Torrent
For a comprehensive list of panels that have been developed or tested through this pipeline, as well as updated and documented associated files, please see the following repository: amplitargets.
- Linux or Mac operating system
- R (and packages within Rscripts)
- simple_pop_stats
- CKMRsim
- Optional tool: SNPlift
- Demo analysis (optional): ms_oyster_panel
If you find this tool useful, please cite the original article that uses the tool:
Sutherland et al. 2024. An amplicon panel for high-throughput and low-cost genotyping of Pacific oyster. G3: Genes|Genomes|Genetics, jkae125, link.
Please also be sure to cite the tools applied within each function.
A. VariantCaller to genepop
B. Population genetic analyses
C. Parentage analysis
D. Database storage
de novo variant calling.
Microhaplotype workflow.
Merge multiple Torrent VariantCaller VCF files.
Convert positions of VCF to a different genome.
Select only hotspot SNPs, excluding novel from TVC VCF.
Looking for panel design? A separate README is available here.
Original approach for technical replicates? here.
Clone this repository and change into the main directory.
git clone https://github.com/bensutherland/amplitools.git
cd amplitools
amplitools is mostly run in R. To initialize a session, open and source 01_scripts/00_initiator.R to initiate functions.
Optional: test out the platform with test data:
load_vc(input_folder = "02_input_data", test_only = TRUE)
This should load 8024 rows and 14 columns, with 2079 unique markers, from two different samples, and will save the output as input.list for inspection.
This section will allow you to convert Torrent VariantCaller output files (TVC) to a standard genepop file for downstream analysis.
An example input file is provided: 02_input_data/test_data.xls
Copy any number of input files in the folder 02_input_data.
Important note on filenames: input filename must be alphanumeric characters only connected with hyphens or underscores (i.e., no spaces). Importantly, multiple sequential underscores without other alphanumeric characters breaking up the underscores will cause the script to fail (e.g., this__will_not_work_as_a_filename.xls).
Essential background on TVC formats:
- 'No Call' = missing data (0000)
- 'Absent' = homozygous reference (0101)
- 'Heterozygous' = heterozygote (0102)
- 'Homozygous' = homozygous alternate (0202)
Where REF or ALT is based on their designations in the hotspot file.
All variantCaller input files must have been generated using the same hotspot file, because this script assumes that all REF or ALT designations as above are the same for all files.
Use the following R function to select only hotspot markers in all input files in 02_input_data, and then convert genotype calls to genepop format to output a multilocus genotype matrix where rows are samples and columns are loci:
proton_to_genepop(neg_control="BLANK")
Function notes:
- sample identifiers will be created as
RunName__Barcode__SampleName. - negative control designation requires the exact matching string
- periods in locus names will be converted to underscores.
- to avoid buggy behaviour, any instance of the string 'pop' in locus names will be converted to '.p.o.p.' to avoid functions that read genepop files into R from seeing this as the separator line between loci names and genotype data (i.e., the POP separator)
- any locus names with a comma in the name will be substituted for an underscore to avoid bugs
Output: tab-delim text file per input file in 02_input_data/prepped_matrices/
To inspect mean marker coverage per experimental or control sample, see output files:
03_results/mean_marker_cov_per_sample_*.txt
From the terminal, finalize the genepop by running the following script for each genotype block text file:
./01_scripts/format_genepop.sh 02_input_data/prepped_matrices/<filename>.txt
The output will be a genepop file for each file output as 02_input_data/prepped_genepops/*.gen
If you have multiple genepops in the amplitools folder, use the following to save them to a list:
# Load the genepops in the input folder
load_genepops(genepop_folder = "02_input_data/prepped_genepops/"
, datatype = "SNP"
, shorten_name = T
)
# Save the list as Rdata to load later in case you clear the workspace
saveRDS(object = my_genepops.list, file = "03_results/my_genepops_list.RDS")
If you are using a single genepop, you can go directly to simple_pop_stats and use the functions there (see below).
Example using simple_pop_stats; to get started:
# clone simple_pop_stats at same level as amplitools
git clone https://github.com/bensutherland/simple_pop_stats.git
cd simple_pop_stats
Copy the genepop(s) created above to simple_pop_stats/02_input_data or analyze elsewhere.
An example is given here.
simple_pop_stats is being developed to work well with amplitools to facilitate population genetic analyses. As follows are some examples of useful functions, but please see the simple_pop_stats repository for a full explanation.
Source simple_pop_stats.
Load the genepop generated by amplitools:
load_genepop(datatype = "SNP")
Simplify your individual names from the amplitools format (i.e., <run>__<barcode>__<indiv>) to simply the indiv ID:
simplify_names(df = obj, format = "amplitools")
Note: must use this function for downstream functions to work.
Create a sample information file to be used for downstream analyses:
generate_popmap(df = obj)
This will produce a form to complete in simple_pop_stats/00_archive/my_data_ind-to-pop.txt
It will automatically produce a row per sample (indiv), and provide empty fields for pop, an alternate identifier (alt.ID), sex, and whether the sample should be ignored.
Manually annotate the above file and save it as my_data_ind-to-pop_annot.txt.
Note: do not include any spaces in field names, and do not start populations with numbers.
Load and population per-sample information based on the completed file above:
annotate_from_popmap(df = obj, popmap.FN = "00_archive/my_data_ind-to-pop_annot.txt")
See simple_pop_stats for other methods, such as filtering, plotting data, analyses, and converting from genepop to rubias format. The below parentage analysis will depend on a rubias file.
To convert the genind to rubias (brief), in R:
# Prepare a stock code file
stock_code.df <- as.data.frame(unique(pop(obj)))
colnames(stock_code.df) <- "collection"
stock_code.df$repunit <- stock_code.df$collection
stock_code.df
# Write it out
write_delim(x = stock_code.df, file = "00_archive/stock_code.txt", delim = "\t", col_names = T)
micro_stock_code.FN <- "00_archive/stock_code.txt"
## All filtered loci: write to rubias for parentage assignment
pop_map.FN <- "00_archive/my_data_ind-to-pop_annot.txt"
genepop_to_rubias_SNP(data = obj, sample_type = "reference"
, custom_format = TRUE, micro_stock_code.FN = micro_stock_code.FN
, pop_map.FN = pop_map.FN
)
# output will be 03_results/rubias_output_SNP.txt
# rename the output and copy it to amplitools/03_results/
The parentage analysis is largely dependent on the R package CKMRsim by Eric C. Anderson. If you use the function ckmr_from_rubias.R, please cite CKMRsim.
Please also see instructions in the associated pipeline: ms_oyster_panel.
Quick method, as an example:
# source amplitools
# Set user variables
input_rubias.FN <- "03_results/rubias_125_ind_583_loc_2024-03-27.txt"
parent_pop <- "VIU_F1"
offspring_pop <- "VIU_F2"
cutoff <- 5
# Prepare run folder
date <- format(Sys.time(), "%Y-%m-%d")
input_rubias_short.FN <- gsub(pattern = "03_results/", replacement = "", x = input_rubias.FN)
input_rubias_short.FN <- gsub(pattern = ".txt", replacement = "", x = input_rubias_short.FN)
run_folder.FN <- paste0("03_results/ckmr_input_", input_rubias_short.FN, "_"
, offspring_pop, "_vs_", parent_pop,"_", date
)
print("Making new result folder...")
print(run_folder.FN)
dir.create(run_folder.FN)
# Run ckmr on the input file
ckmr_from_rubias(input.FN = input_rubias.FN
, parent_pop = parent_pop
, offspring_pop = offspring_pop
, cutoff = cutoff
, output.dir = run_folder.FN
)
It is now possible to graph the relatives based on the output of CKMR-sim, using the following function, with inputs being the paired relatives and associated log likelihood values (.txt format):
graph_relatives(input.FN = "03_results/parent_fs_broodstock_pw_logl_5.txt"
, logl_cutoff = 5 # Additional logl cutoff as needed
, drop_string = "G00" # String constant to remove from sample names
, directed = FALSE # Should the plotted graph be directed?
, plot_width = 5 # Plotting width
, plot_height = 5 # Plotting height
)
Other inputs may include offspring sibship, or parent-offspring relationships.
This section provides suggestions to maintain amplicon sequence data from variantCaller outputs for long-term storage.
We suggest the following approach:
- keep all raw data stored separately from amplitools;
- copy in a specific raw data runs to a fresh instance of the main repo;
- save an output file in standard genetics format (e.g., rubias) as a 'curated dataset'. This will be date-stamped, and then should not be altered. If it needs to be altered, return to the start and re-run.
Here are several general tips on how to manage this:
As an example considering Pacific oyster, on the server/ computer, in a safe and permanent location, do the following:
mkdir 00_raw_data
cd 00_raw_data
mkdir cgig
cd cgig
mkdir cgig_v.1.0
Once the structure is built, copy raw run files (.xls) into the above folder.
- within the folder, store a single spreadsheet explaining per chip the chip name, date run, and projects included. This will be used to track what is in each raw data file, so specific chips can be then copied into an active amplitools repo.
- note: it may be best to make the above read-only without sudo priviledges, to be sure that files do not get altered
- also include an md5 alongside each raw run file, with the same name but .md5 appended to the end of the filename.
- also include a sample interp file alongside each chip (see below).
e.g., R_<your_chip_name>_sample_interp.txt
Include the following columns so this file can also be used as a popmap:
| indiv | pop | alt.ID | sex | ignore |
|---|---|---|---|---|
| F2-03 | VIU_F2 | NA | NA | NA |
| BR33 | VIU_F1 | NA | NA | NA |
| ... | ... | ... | ... | ... |
Note: the ignore column will eventually be used so that one can set amplitools to ignore the sample, for example if it is decided that the sample should not be included in building downstream datasets, or perhaps it is the poorer of two technical replicates and should be dropped. This is currently not implemented, but is planned in a future release.
Importantly, when a new hotspot is used to generate the output, the files will no longer be directly comparable, and so a new folder should be created (e.g., cgig_v.1.1) for clarity.