---
title: "12.1-snp-calling-BS-Snper-parents"
author: "Kathleen Durkin"
date: "`r Sys.Date()`"
always_allow_html: true
output: 
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: true 
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
---

To mask C/T SNPs during methylation calling work and to do any relatedness/structure work with the WGBS data, I need to first call SNPs using a methylation-aware caller (in this case, BS-Snper). I believe Venkataraman et al. 2024 already did this for the parent data, but the only files I could find on the `ceabigr` repo were for C/T SNPs -- I also want other SNPs for the relatedness work.

To get those, and to simply ensure consistent parameters/tool versions/etc. across the offspring and parent SNPs, I'll be re-calling parent SNPs here.

I'll be using `BS-Snper` for WGBS variant calling, and referencing a past example from Yaamani (accessed through the [Roberts Lab Handbook](https://robertslab.github.io/resources/bio_DNA-methylation/#bs-snper))


# Load libraries

```{r}
library(dplyr)
library(tidyr)
library(pheatmap)
```


# Download all sample BAM files 

# Download methylation calls

The  `bismark` .deduplicated.sorted.bam coverage files for the ceabigr parent gametes stored at: https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/120321-cvBS/

Importantly, not all of these individuals were included in the crosses that generated our study's offspring individuals. I need to ID the parents of all of our individuals, get their original sample names, and then download their respective files. 

```{r download-parent-files, engine='bash', eval=FALSE}

# For all offspring, grab their two parent IDs from filenames
parents=$(ls ../output/02.10-bismark-deduplication/*dedup.sorted.bam \
| xargs -n1 basename \
| sed -E 's/^([A-Za-z0-9]+)-([A-Za-z0-9]+)-.*/\1\n\2/' \
| sort -u)

# from parent IDs, identify sex (M/F) and treatment (C/E)
# also, from metadata file, use Sample ID to get original Parent ID (used in parent data file names)
: > "../data/bismark-methyl-extraction-parents/parent_lookup.tsv"
printf 'Parent.ID\ttreatment\tsex\tSample.ID\n' >> "../data/bismark-methyl-extraction-parents/parent_lookup.tsv"
sample_ids=()
while read -r p; do
  case ${p:0:1} in C) trt=Control;; E) trt=Exposed;; *) trt=NA;; esac
  case ${p:1:1} in F) sex=Female;;  M) sex=Male;;    *) sex=NA;; esac
  sid=$(awk -F, -v pid="$p" 'NR>1 && $6==pid {print $1; exit}' "../data/adult-meta.csv")
  if [ -z "$sid" ]; then
    printf '%s\t%s\t%s\t<NO MATCH>\n' "$p" "$trt" "$sex" >> "../data/bismark-methyl-extraction-parents/parent_lookup.tsv"
    echo "WARNING: no metadata match for parent $p , skipping" >&2
  else
    printf '%s\t%s\t%s\t%s\n' "$p" "$trt" "$sex" "$sid" >> "../data/bismark-methyl-extraction-parents/parent_lookup.tsv"
    sample_ids+=("$sid")
  fi
done <<< "$parents"

# now that we have list of Parent IDs, download corresponding .cov.gz files from storage
alt=$(IFS='|'; echo "${sample_ids[*]}")
wget -r -np -nd -nc \
     --accept-regex "/(${alt})[^/]*\.deduplicated\.sorted\.bam$" \
     -P "../data/bismark-methyl-extraction-parents" "https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/120321-cvBS/"

```

```{r check-download, engine='bash', eval=FALSE}
cd ../data/bismark-methyl-extraction-parents

echo "How many bam files were downloaded?"
ls *deduplicated.sorted.bam | wc -l

echo ""

echo "Check file checksums:"
grep 'deduplicated.sorted.bam' checksums.md5 | md5sum -c -
```


# Merge samples

First, to ID SNPs *accross* all samples, for use in removing C/T transitions during methyl analysis, I first need to merge all sample BAMs into a single BAM file, which will be used as input to `BS-Snper`.

merge to single BAM using `samtools`
```{r, engine='bash', eval=FALSE}
samtools merge \
--threads 20 \
-o ../output/12.1-snp-calling-BS-Snper-parents/merged_offspring_deduplicated.sorted.bam \
../data/bismark-methyl-extraction-parents/*deduplicated.sorted.bam \
```

check
```{r, engine='bash', eval=FALSE}
samtools view \
--threads 20 \
../output/12.1-snp-calling-BS-Snper-parents/merged_offspring_threads_deduplicated.sorted.bam \
| head
```

create index file (in case I want to view in IGV)
```{r, engine='bash', eval=FALSE}
samtools index \
-b ../output/12.1-snp-calling-BS-Snper-parents/merged_offspring_deduplicated.sorted.bam 
```


# Identify SNPs across samples

Options for the script are found here and below. Input and output files should be absolute paths.

--fa: Reference genome file in fasta format
--input: Input bam file (I'm using deduplicated sorted bams)
--output: Temporary file storing SNP candidates
--methcg: CpG methylation information
--methchg: CHG methylation information
--methchh: CHH methylation information
--minhetfreq: Threshold of frequency for calling heterozygous SNP
--minhomfreq: Threshold of frequency for calling homozygous SNP
--minquali: Threshold of base quality
--mincover: Threshold of minimum depth of covered reads
--maxcover: Threshold of maximum depth of covered reads
--minread2: Minimum mutation reads number
--errorate: Minimum mutation rate
--mapvalue: Minimum read mapping value
SNP.out: Final SNP result file
ERR.log: Log file

```{r, engine='bash', eval=FALSE}
#Defaults: --minhetfreq 0.1 --minhomfreq 0.85 --minquali 15 --maxcover 1000 --minread2 2 --errorate 0.02 --mapvalue 20
#Saving output to gannet
/home/shared/BS-Snper-master/BS-Snper.pl \
--fa ../data/Cvirginica_v300/Cvirginica_v300.fa \
--input ../output/12.1-snp-calling-BS-Snper-parents/merged_offspring_deduplicated.sorted.bam \
--output ../output/12.1-snp-calling-BS-Snper-parents/SNP-candidates.txt \
--methcg ../output/12.1-snp-calling-BS-Snper-parents/CpG-meth-info.tab \
--methchg ../output/12.1-snp-calling-BS-Snper-parents/CHG-meth-info.tab \
--methchh ../output/12.1-snp-calling-BS-Snper-parents/CHH-meth-info.tab \
--mincover 5 \
> ../output/12.1-snp-calling-BS-Snper-parents/SNP-results.vcf 2> ../output/12.1-snp-calling-BS-Snper-parents/merged.ERR.log

```







# Identify SNPs per-sample

To call per-sample variants (to ultimately generate a relatedness matrix), I need to run `BS-Snper` separately one each sample BAM. Since there are >30 samples, want to parallelize this to save time.
```{r, engine='bash', eval=FALSE}
mkdir -p ../output/12.1-snp-calling-BS-Snper-parents/per-sample

# function to run BS-Snper once
run_bssnper() {
  bam=$1
  s=$(basename ${bam} deduplicated.sorted.bam)
  ~/BS-Snper-master/BS-Snper.pl \
    ${bam} \
    --fa ../data/Cvirginica_v300/Cvirginica_v300.fa \
    --output ../output/12.1-snp-calling-BS-Snper-parents/per-sample/${s}.candidates.txt \
    --methcg ../output/12.1-snp-calling-BS-Snper-parents/per-sample/${s}.CpG.tab \
    --methchg ../output/12.1-snp-calling-BS-Snper-parents/per-sample/${s}.CHG.tab \
    --methchh ../output/12.1-snp-calling-BS-Snper-parents/per-sample/${s}.CHH.tab \
    --mincover 5 \
    > ../output/12.1-snp-calling-BS-Snper-parents/per-sample/${s}.vcf 2> ../output/12.1-snp-calling-BS-Snper-parents/per-sample/${s}.err.log \
  && rm ../output/12.1-snp-calling-BS-Snper-parents/per-sample/${s}.candidates.txt ../output/12.1-snp-calling-BS-Snper-parents/per-sample/${s}.C*.tab
}
export -f run_bssnper

# list out all of the sample BAM files, then feed into function run_bssnper using `parallel` to parallelize
ls ../data/bismark-methyl-extraction-parents/*deduplicated.sorted.bam | parallel -j 20 --joblog ../output/12.1-snp-calling-BS-Snper-parents/per-sample/parallel.log run_bssnper {}

```


# Isolate C->T SNPs

I'll need to mask C->T SNPs before methylation analysis, so save file of just C/T SNPs for each sample.

Note, however, that this will actually mean masking both C->T SNPs *and* G->A SNPs, to account for strandedness (a C->T SNP on the - strand would look like a G-> SNP on the + strand)

```{r, engine='bash'}
cd ../output/12.1-snp-calling-BS-Snper-parents/per-sample

for file in *.vcf
do
  # grep -v "^#" skips the header lines
  # the awk selects only lines with either a C/T or G/A in columns 4/5
  grep -v "^#" ${file} \
  | awk -F'\t' '($4=="C" && $5=="T") || ($4=="G" && $5=="A"))' \
  > ${file}_CT_GA_SNPs.tab
  
  # check files
  head ${file}_CT_GA_SNPs.tab
  wc -l ${file}_CT_GA_SNPs.tab
done

```

# Merge vcf files

Now that I have a variant file for each sample, combine into single file for conversion to relatedness matrix

File should end up as a multi-sample .vcf.gz file, which has one row per SNP position, one genotype column per sample.
```{r, engine='bash', eval=FALSE}
cd ../output/12.1-snp-calling-BS-Snper-parents/per-sample

# For each sample's output .vcf file, pull the sample name and rewrite file header
# to include unique sample ID and chromosome names/lengths (required for BS-Snper)
# Also compress to .vcf.gz file type, which is required to merge all files using `bcftools` merge
for vcf in *.vcf; do
  s=${vcf%.vcf}
  echo ${s} > ${s}.name
  /home/shared/bcftools-1.14/bcftools reheader -s ${s}.name --fai ../../../data/Cvirginica_v300.fa.fai ${vcf} \
    | /home/shared/bcftools-1.14/bcftools view -f PASS -Oz -o ${s}.vcf.gz
  /home/shared/bcftools-1.14/bcftools index ${s}.vcf.gz
done

# merge all files using `bcftools` merge
ls *.vcf.gz > vcf_list.txt
/home/shared/bcftools-1.14/bcftools merge -l vcf_list.txt -Oz -o merged.vcf.gz
```


Filter to keep: 

- only biallelic sites
- only single base substitutions (SNPs)
- only sites with genotypes in at least half of the samples
- only sites where the rarer allele is at a frequency of at least 5% (rare variants not super useful for relatedness estimaes, add noise)
```{r, engine='bash', eval=FALSE}
cd ../output/12.1-snp-calling-BS-Snper-parents/per-sample

/home/shared/bcftools-1.14/bcftools view -m2 -M2 -v snps merged.vcf.gz \
  | /home/shared/bcftools-1.14/bcftools view -i 'F_MISSING < 0.5 && MAF > 0.05' -Oz -o merged.filtered.vcf.gz
```

Use `plink` to build the relatedness matrix (`MACAU` reccommends using `GEMMA`, but we already have `plink` installed on Raven at `/home/shared/plink_linux_x86_64_20230116/plink`)


`--allow-extra-chr` : `plink` expects human-style chromosomes, need this option to permit the oyster chromosomes
```{r, engine='bash', eval=FALSE}
/home/shared/plink_linux_x86_64_20230116/plink --vcf merged.filtered.vcf.gz \
  --allow-extra-chr --double-id \
  --set-missing-var-ids @:# \
  --make-rel square \
  --out relatedness
```


Build a heatmap showing the sample relatedness
```{r}
rel_dir <- "../output/12.1-snp-calling-BS-Snper-parents/per-sample"

# read in relatedness matrix
rel <- as.matrix(read.table(file.path(rel_dir, "relatedness.rel")))

# read sample ids from the helper doc
ids <- read.table(file.path(rel_dir, "relatedness.rel.id"))$V2
dimnames(rel) <- list(ids, ids)

# Parse annotations from sample names
# (e.g., "CF01-CM01-Zygote" used to pull Control and Zygote)
parts <- do.call(rbind, strsplit(ids, "-"))
annot <- data.frame(
  Stage     = parts[, 3],
  Treatment = substr(parts[, 1], 1, 1),   # C or E
  row.names = ids
)

# Build heatmap
pheatmap(rel,
         annotation_row = annot,
         annotation_col = annot,
         color = colorRampPalette(c("#2166ac", "white", "#b2182b"))(100),
         breaks = seq(-max(abs(rel)), max(abs(rel)), length.out = 101),
         fontsize_row = 7, fontsize_col = 7,
         main = "Genomic relatedness (PLINK --make-rel)")
```

Looks good! Half-sibs are grouped together, as expected, and grouping is associated with neither treatment nor stage, which is also great! Only confusing thing is why `CF08-CM05-Larvae` and `EF07-EM01-Zygote` are so closely related...

It could just be because one of those samples has extremely low coverage, and they share a relatively high proportion of SNPs by chance (`CF08-CM05-Larvae` only has 6.2 million C's, compared to >200 million in most samples)