HybSeq-SNP-Extraction

Software

HybPiper (produces supercontigs): https://github.com/mossmatters/HybPiper
GATK4: https://github.com/broadinstitute/gatk/releases
Plink: https://www.cog-genomics.org/plink/1.9/
Samtools: https://github.com/samtools/samtools
BWA: https://github.com/lh3/bwa

Additional Software For Analyses

Haplonerate: https://github.com/mossmatters/phyloscripts/tree/master/haplonerate
WhatsHap: http://whatshap.readthedocs.io
BioPython package: https://biopython.org/
Python 3.0: https://www.python.org/download/releases/3.0/
GNU Parallel: https://www.gnu.org/software/parallel/
bcftools: https://samtools.github.io/bcftools/

Prerequisites

From HybPiper, first run the scripts reads_first.py followed by intronerate.py to produce supercontigs for each recovered gene. These supercontigs will be used as a "reference sequence" for the sample. If you have several different species, you will need to run the scripts for each individually.

Concatenate all supercontigs into one single (reference) file:
prefix/*/prefix/sequences/intron/*_supercontig.fasta > prefix.supercontigs.fasta

Note: "prefix.supercontigs.fasta" will be used as an input on command line along with samplename

Workflow

variantcall.sh

This script will:

1. Uses bwa mem to map paired-end reads to supercontigs
2. Replaces read groups for mapped and unmapped bam files
3. Removes duplicate reads
4. Identifies variant sites using gatk HaplotypeCaller (NOTE: GVCF is produced)
5. Removes intermediate BAM files

Command line: bash variantcall.sh prefix.supercontigs.fasta samplename

Output: Contains many intermediate BAM files and GVCF file

GenotypesToPCA.sh

This script will:

1. Create samples.list from GVCF files (Use samples.list as variant in step 2)
2. Combine GVCF files into a cohort and genotype
3. Filter SNP's to remove indels using hard filter
   "QD < 5.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"
4. Select SNP variants from hard filtering parameters

Command line: bash GenotypesToPCA.sh prefix.supercontigs.fasta species

Potentially useful outputs: "$prefix".SNPall.vcf (contains all SNPs and indels), "$prefix".snp.filtered.vcf (only SNPs, removes indels), "$prefix".snp.filtered.nocall.vcf (ONLY SNPs that pass a hard filter)

plink_stats.sh

Additional dependencies:

bcftools: https://samtools.github.io/bcftools/
Plink: https://zzz.bwh.harvard.edu/plink/download.shtml

This script will:

1. Set ID name for each SNP (for filtering)
2. Filter SNPs that didn't pass the filter or have missing data
3. Generate eigenvalues and loadings for PCA axes (default set to 20)
4. Generate basic statistics (heterozygosity, inbreeding coefficient, allele frequencies)

Command line: bash plink_stats.sh "$prefix"

Outputs: Statistic files "$prefix"_pruned.bed, "$prefix"_pruned.bim, "$prefix"_pruned.fam, "$prefix"_pruned.frq, "$prefix"_pruned.het, "$prefix"_pruned.ibc, "$prefix"_pruned.map, "$prefix"_pruned.nosex, "$prefix"_pruned.ped, "$prefix"_pruned.log

extract_phase_subgenomes.sh

NOTE: DO NOT RUN HAPLOTYPECALLER IN GVCF MODE

This workflow has been modified from Kates et al paper for allodiploid genomes
Link to Kates et al paper: https://pubmed.ncbi.nlm.nih.gov/29729187/

Additional Dependencies:

Haplonerate: https://github.com/mossmatters/phyloscripts/tree/master/haplonerate
WhatsHap: http://whatshap.readthedocs.io
BioPython package: https://biopython.org/
Python 3.0: https://www.python.org/download/releases/3.0/
GNU Parallel: https://www.gnu.org/software/parallel/

This script will:

1. Replace the FASTA headers
2. Run WhatsHap to generate phased VCF using pipe characters
3. Extract two FASTA sequences for each gene, corresponding to two alleles

Command line: bash extract_phase_subgenomes.sh "$prefix"
Output: One file containing two FASTA format subgenome sequences per individual per gene to be used in phylogenetic analysis