Bioinformatics manual for population epigenomics combining whole-genome and target genome sequencing

java -Xmx4G -jar trimmomatic.jar PE -threads 12 file_R1.fastq.gz file_R2.fastq.gz 
file_trimmed_1.fastq.gz  file_unpaired_1.fastq.gz file_trimmed_2.fastq.gz
file_unpaired_2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:35

bwa mem genome.fa file_trimmed_1.fastq.gz file_trimmed_2.fastq.gz -t 12 -M > file.sam 

samtools view -Sb file_trimmed.sam > file_trimmed.bam 
samtools sort file_trimmed.bam  -o file_trimmed_sorted.bam 
samtools flagstat file_trimmed_sorted.bam > file_flagstats.txt 
samtools stats file_trimmed_sorted.bam > file_stats.txt 

interleave-reads.py file_R1.fastq file_R2.fastq -o file_interleave_R1_R2.fastq

normalize-by-median.py -k 20 -C 30 -N 4 -x 4e9 file_interleave_R1_R2.fastq  -o file_normalize_by_median_R1_R2.fastq

extract-paired-reads.py file_normalize_by_median_R1_R2.fastq -f --output-paired file_diginorm_paired --output-single file_diginorm_single

java -Xmx16g -jar picard.jar MarkDuplicates I=file_trimmed_sorted.bam O=file_trimmed_sorted_rmdup.bam CREATE_INDEX=true REMOVE_DUPLICATES=true M=file_output.metrics

## HaplotypeCaller
gatk --java-options "-Xmx16G" HaplotypeCaller -R genome.fa -I file_trimmed_sorted_rmdup.bam -ERC GVCF -O file_trimmed_sorted_rmdup.g.vcf
## GenomicsDBImport
gatk --java-options "-Xmx96G -Xms96G" GenomicsDBImport -V file1_trimmed_sorted_rmdup.g.vcf -V file2_trimmed_sorted_rmdup.g.vcf --genomicsdb-workspace-path my_database -L list_Chr+scaff.list --batch-size 50 -ip 500
## GenotypeGVCFs
gatk GenotypeGVCFs -R genome.fa -V gendb://my_database  -new-qual true -O all_trimmed_sorted_rmdup_gVCF_GATK.snps.indels.vcf

#HaplotypeCaller
GATK -R haplome_v2.3.fa -T HaplotypeCaller -nct 20 -I sample1_trimmed_vs_haploV23.bam -I sample2_trimmed_vs_haploV23.bam -I sample3_trimmed_vs_haploV23.bam -I sample4_trimmed_vs_haploV23.bam -I
 sample5_trimmed_vs_haploV23.bam -I sample6_trimmed_vs_haploV23.bam -I sample7_trimmed_vs_haploV23.bam -I sample8_trimmed_vs_haploV23.bam -I sample9_trimmed_vs_haploV23.bam -I sample9_trimmed_vs_haploV23.bam 
--emitRefConfidence GVCF -o gatk_nct20_slurm_1node-c20_snps.vcf

#GenotypeGVCFs
GATK -T GenotypeGVCFs -R haplome_v2.3.fa --variant sample1.vcf --variant sample2.vcf --variant sample3.vcf --variant sample4.
vcf --variant sample5.vcf --variant sample6.vcf --variant sample7.vcf --variant sample8.vcf --variant sample9.vcf --variant sample10.vcf -o gatk_all10samples_SNPs.vcf

samtools mpileup -uf genome.fa mapping_file_sort_without_duplicate.bam | bcftools call -mv -Oz > file_bcftools_noduplicate.vcf.gz

freebayes -f genome.fa all_samples.bam > freebayes_all_samples.vcf

vcftools --vcf all_tool.snps.indels.vcf --out all_filtered_tool.vcf --remove-indels --max-alleles 2 --min-alleles 2 --minQ 30--recode --recode-INFO-all

bcftools index sample1_diginorm_gatk3.8_depth20_maf30.vcf.gz
bcftools index sample1_diginorm_FreeBayes_depth20_maf30.vcf.gz
bcftools index sample1_samtools_depth20_maf30.vcf.gz

bcftools isec -n +3 sample1_diginorm_gatk3.8_depth20_maf30.vcf.gz sample1_diginorm_FreeBayes_depth20_maf30.vcf.gz sample1_samtools_depth20_maf30.vcf.gz -O v -o common_SNPs_sample1_GATK_FreeBayes_samtools_depth20_maf30_bcftools.txt

methratio.py   --ref ref_genome.fa --zero-meth TRUE --trim-fillin 2 --combine-CpG --min-depth 8 --context all  bsmap_sample*.sam

meth.CpG <- unite(CpG, destrand = TRUE, min.per.group = 7L)
meth.CHG <- unite(CHG, destrand = FALSE, min.per.group = 7L)
meth.CHH <- unite(CHH, destrand = FALSE, min.per.group = 7L)

SNPdat <- read.delim("SNP_file.txt", header = F)
     
#with SNP_file.txt:
#    ScaffoldID     position     allele1     allele2

SNPdat$Scaff_Pos <- paste(SNPdat$Scaff, SNPdat$Pos, sep="_")
SNPdat$SNP <- paste(SNPdat$Ref, SNPdat$Alt, sep ="/")
MethPos2 <- paste(meth.CpG2$chr, meth.CpG2$start, sep = "_")
MethPosMatchSNP2 <-which(MethPos2 %in% SNPdat$Scaff_Pos)
SNPMeth2 <- subset(SNPdat, Scaff_Pos %in% MethPos2[MethPosMatchSNP2])
SNPMethOk <- subset(SNPMeth2, SNP == "C/T")
CpG.posOK2 <- select(meth.CpG2, which (!MethPos2 %in% SNPMethOk$Scaff_Pos)) 

for (i in 1:19) {
cov <- getData(meth.CHG.filtind.filtSNP.filtCov)[,colnames(meth.CHG.filtind.filtSNP.filtCov) == paste0("coverage", i)]
cov_filt <- sort(c(which(cov < 7), which(is.na(cov)))) meth.CHG.filtind.filtSNP.filtCov[cov_filt, colnames(meth.CHG.filtind.filtSNP.filtCov) == paste0("numCs", i)] <- NA meth.CHG.filtind.filtSNP.filtCov[cov_filt, colnames(meth.CHG.filtind.filtSNP.filtCov) == paste0("numTs", i)] <- NA
     rm(cov, cov_filt)
}

meth.CpG.window <-tileMethylCounts(meth.CpG.filtind.filtSNP.filtTE.filtCov.filtNA,win.size = 1000, step.size = 250)
meth.CHG.window <-tileMethylCounts(meth.CHG.filtind.filtSNP.filtTE.filtCov.filtNA,win.size = 1000, step.size = 250)
meth.CHH.window <-tileMethylCounts(meth.CHH.filtind.filtSNP.filtTE.filtCov.filtNA,win.size = 1000, step.size = 250)

#Identification of windows to remove
percmeth.CpG.window.sd <- rowSds(percmeth.CpG.window, na.rm = TRUE)
sum(percmeth.CpG.window.sd == 0)

# Removal of windows showing the less variable levels of methylation
percmeth.CpG.window <-percmeth.CpG.window[which(percmeth.CpG.window.sd != 0), ] dim(percmeth.CpG.window)

#Identification of the windows associated with the most variable methylation levels
percmeth.CpG.window.sd <- rowSds(percmeth.CpG.window, na.rm = TRUE)
layout(matrix(c(rep(1, 2), 2), nrow = 1))
hist(percmeth.CpG.window.sd, col = "grey", main = "")
bp <- boxplot(percmeth.CpG.window.sd, col = "grey")
length(bp$out)
bp$stats

dag_window_size=1000
dag_step=250

load("meth.CHG.filtind.filtSNP.filtTE.filtCov.filtNA.Rdata")
y<-x[,c("chr","start","end","strand")]

for (i in 1:length(colnames(x)[colnames(x) %like% "coverage"])){  # To recover the C/coverage values
  j=5+3*(i-1)
  print(paste0(j," ",j+1))    
  y[,paste0("in",i)]<-x[,j+1]/x[,j]
}
yy<-x[,c("chr","start","end","strand")]
rm(x)

z<-rowSds(as.matrix(y[,5:ncol(y)]),na.rm=TRUE) # Calculate row standard deviations
yy$STDEV<-z
rm(z)
y<-yy
rm(yy)


# Do last adaptations and launch
dag_window=dag_window_size/dag_step
colnames(y)<-c("CHR","START","END","STRAND","STDEV")
y$MEAN<-(y$START+y$END)/2
y$CHR<-gsub("Chr0","Chr",y$CHR,perl=TRUE)
y$WINDOW<-floor(y$MEAN/dag_step)+1


stdev_counts = data.table(
   CHR = character(),
   WIN = numeric(),
   POS = numeric(),
   STDEV = numeric()
)

count=0
for (i in unique(y[y$CHR %like% "Chr" | y$CHR %like% "scaffold",]$CHR)){
  window_size=dag_window_size
  step=dag_step
  #i<-paste0("Chr",i)
  z<-y[y$CHR==i,]
  min=0
  max=max(z$WINDOW)
  #print(paste(i,min,max,min(z$MEAN),max(z$MEAN)))
  count=count+1
  print(paste(i,min,max,min(z$MEAN),max(z$MEAN),count,length(unique(y[y$CHR %like% "Chr" | y$CHR %like% "scaffold",]$CHR))))
  zz<-data.frame(matrix(ncol=2,nrow=max*step))
  colnames(zz)<-c("MEAN","STDEV")
  zz$MEAN<-rownames(zz)
  
  zz[zz$MEAN %in% z$MEAN,]$STDEV<-z[z$MEAN %in% zz$MEAN,]$STDEV

 # Sliding window
  total <- nrow(zz)
  if (max(z$MEAN)<window_size){ # Adapted to avoid problems with scaffolds smaller than window_size
    spots <- 1
  }
  else {
    spots <- seq(from=1, to=(total-window_size), by=step)
  }
  
  if (spots[length(spots)]<=total-window_size){spots<-c(spots,(spots[length(spots)]+step))} # Adapted to recover the last bits inside smaller window
  result <- vector(length = length(spots))
  for(j in 1:length(spots)){
    if (j%%50000==0){print(paste(j,length(spots)))}
    if ((spots[j]+window_size)>=total){window_size=(total-spots[j])} # Adapted to recover the last bits inside last smaller window 
    result[j] <- mean(zz[spots[j]:(spots[j]+window_size-1),"STDEV"],na.rm=TRUE)
  }
  
  stdev_counts<-rbind(stdev_counts,data.frame(CHR=i,WIN=1:length(spots),POS=spots,STDEV=result))
}

x<-stdev_counts
write.table(x,file=paste0(save_file_name))

#$Id$

###run with python get_threshold_over_all_windows_calc1.py OUTPUT_FILE_from_calc1_get_mean_and_stdv_for_each_window.py > threshold_calc1.txt


import os
import re
import string
import sys
import glob
import numpy

file1 = sys.argv[1]
file1_stream = open(file1)
list_of_means = []

for line1 in file1_stream.readlines():
	if (line1.count('start') == 0):
		line1 = line1.replace('\n','')
		splitted_line1 = line1.split('\t')
		scaffold = splitted_line1[0]
		start = splitted_line1[1]
		end = splitted_line1[2]
		
		mean = splitted_line1[13]
		mean = float(mean)
		list_of_means.append(mean)

list_of_means.sort()
nbre_de_means = len(list_of_means)
##XXX corresponds to the first half of the dataset 
##YYY corresponds to the second half of the dataset
Q1 = numpy.median(list_of_means[:XXX])
Q3 = numpy.median(list_of_means[YYY:])

##for CHH context, hreshold = (Q3 + 3*(Q3- Q1))
threshold = (Q3 + 1.5*(Q3- Q1))
threshold = round(threshold,5)

print 'threshold = ',threshold

#$Id$

###run with python get_threshold_stdv_over_all_windows_calc2.py OUTPUT_FILE_from_get_stdv_between_individuals_for_each_window_calc2.py > threshold_calc2.txt


import os
import re
import string
import sys
import glob
import numpy

file1 = sys.argv[1]
file1_stream = open(file1)
list_of_stdv = []

for line1 in file1_stream.readlines():
	if (line1.count('start') == 0):
		line1 = line1.replace('\n','')
		splitted_line1 = line1.split('\t')
		scaffold = splitted_line1[0]
		start = splitted_line1[1]
		end = splitted_line1[2]
		
		stdv = splitted_line1[4]
		stdv = float(stdv)
		list_of_stdv.append(stdv)

list_of_stdv.sort()
nbre_de_stdv = len(list_of_stdv)
##XXX corresponds to the first half of the dataset 
##YYY corresponds to the second half of the dataset
Q1 = numpy.median(list_of_stdv[:XXX])
Q3 = numpy.median(list_of_stdv[YYY:])

##for CHH context, hreshold = (Q3 + 3*(Q3- Q1))
threshold = (Q3 + 1.5*(Q3- Q1))
threshold = round(threshold,5)

print 'threshold = ',threshold

trim_galore input_R1.fastq.gz input_R2.fastq.gz --paired ADAPTER1 -a2 ADAPTER2 -o output_directory --gzip -j {threads} 

fastqc trimmed_reads.fq.gz -o fastQC_output_directory -t {threads}

bsmapz -a fileR1.fq.gz -b fileR2.fq.gz -o {output.out} -d mapped_file.bam -d ref_genome.fa -p threads
 

samtools stats sample_bsmapz_sorted.bam -r ref_genome.fa -@ {threads} > sample.statics
samtools fixmate -@ {threads} -O BAM -m  sample_bsmapz_sorted.bam  sample_fixmate.bam 
samtools sort -@ {threads} -O BAM sample_fixmate.bam  -o sample_fixmate_sort.bam 
samtools markdup -r  ref_genome.fa -@ {threads} -s -f sample.statics sample_fixmate_sort.bam sample_fixmate_sort_temp.bam

python methratio.py sample.dedup.bam -o meth_sample.txt -d ref_genome.fa -N {threads} -I

SeqCapBis_CHG = methRead(location = path_to_the_files, sample.id = sample.ids, assembly = "quercus", mincov = 10, context = "CHG", treatment = rep(0,10))

#!/bin/bash 
# Splitting context:

usage()
{
cat << EOF
usage: $0 <options>
splitting context.
OPTION:
  -h      show this Help message.
  -o      Output.
  -i      Input.
EOF
}

# Get options
while getopts "ho:i:" OPTION
do
  case $OPTION in
    h)  usage; exit 1;;
    o)  output=$OPTARG;;
    i)  input=$OPTARG;;
    ?)  usage; exit;;
  esac
done

# Check that all options were passed
if [[ -z $output ]] || [[ -z $input ]] 
then
  printf "\n=========================\n ERROR: missing options\n=========================\n\n"
  usage
  exit 1
fi

#in_file = snakemake.input["isoforms"]
#out_file = snakemake.output["plot"]

# Fail on the first error
set -e

######################

file=$(echo $output|rev|cut -d "/" -f 1 |rev)
path=$(echo $output|rev|cut -d "/" -f 2- |rev)

for context in "CHH" "CG" "CHG"; do 
   
   awk "NR<=1 || \$4~/$context/" $input  > $path/$context-$file ;
done

# Read methylation using methylkit function methRead
myobj <- methRead(location = files, sample.id = sample_id, assembly = "populus tricharpa v3.1", mincov = 1, context = context,treatment = rep(0, length(files)), pipeline = list(fraction=TRUE, chr.col=1, start.col=2, end.col=2, coverage.col=6, strand.col=3, freqC.col=5 ))  

# Concatenate all samples tables into one unique table
finalFrame <- mergeMethylkitOutput(myobj)

#Write the final table as a csv2 file
write.csv2(finalFrame,file = table,)

# head(myobj)

# plots for statistcs and coverage simple :
pdf(file = XXX)
getMethylationStats(myobj[[1]],plot=TRUE,both.strands=FALSE)
getCoverageStats(myobj[[1]],plot=TRUE,both.strands=FALSE)
dev.off()

#Trim the paired sequences
trim_galore -q 30 --paired -o paired_1.fastq  paired_2.fastq 

#Index Genome
bwa index genome_ref.fa

#Align
bwa mem -Y genome_ref.fa  paired_trimmed_1.fastq  paired_trimmed_2.fastq > whole.sam

#From SAM2BAM
samtools view -S -b whole.sam -o whole.bam

#Extract Unmapped reads

#An unmapped read whose mate is mapped.
samtools view -u  -f 4 -F264 whole.bam  > tmps1.bam

#Both reads of the pair are unmapped
samtools view -u -f 12 -F 256 whole.bam > tmps2.bam

#merge
samtools merge unmapped.bam tmps1.bam tmps2.bam

#Extract the reads in FASTQ format (paired)
bamToFastq -bam unmapped.bam -fq1 unmapped_reads1.fastq -fq2 unmapped_reads2.fastq

WD="path/to/working/_directory"
PREFIX="prefix_you_want"

##For each samples
python teflon_prep_custom.py -wd ${WD}reference -g genome_ref -l path/to/TE_LIBRARY -p ${PREFIX}

bwa index ${WD}reference/${PREFIX}.prep_MP/${PREFIX}.mappingRef.fa

bwa mem -Y ${WD}reference/${PREFIX}.prep_MP/${PREFIX}.mappingRef.fa ${READS1} ${READS2} > ${WD}reference/${PREFIX}.sam

samtools view -Sb ${WD}reference/${PREFIX}.sam | samtools sort -o ${WD}reference/${PREFIX}.sorted.bam

samtools index ${WD}reference/${PREFIX}.sorted.bam

#Run Teflon
#For each samples
python teflon.v0.4.py -wd ${WD} -d ${WD}reference/${PREFIX}.prep_TF/ -s path/to/samples -i unique_ID -l1 family -l2 class

#Teflon collapse
##Only once
python teflon_collapse.py -wd ${WD} -d ${WD}reference/${PREFIX}.prep_TF/ -s path/to/samples -n1 minimum_reads_to_support_TE_in_one_sample -n2 minimum_reads_to_support_TE_in_all_samples

#Teflon Count
#For each samples
python teflon_count.py -wd ${WD} -d ${WD}reference/${PREFIX}.prep_TF/ -s path/to/samples -i unique_ID

#Teflon genotype
##Only once
python teflon_genotype.py -wd ${WD} -d ${WD}reference/${PREFIX}.prep_TF/ -s path/to/samples -dt pooled

bismark_genome_preparation --verbose genome_ref.fa

bismark --genome genome_ref.fa paired_trimmed_1.fastq  paired_trimmed_2.fastq  --un

idxEpiTEome.pl -l 100 -gff genome_ref.gff -t /path/to/TE_LIBRARY -fasta genome_ref.fa

epiTEome.pl -gff genome_ref.gff -ref genome_ref.epiTEome.masked.fasta -un unmapped_reads.fastq -t /path/to/TE_LIBRARY