################################################
## Exercice 0 : préparation de l'environnement de travail pour le TP
################################################
## Se positionner dans le répertoire de travail accessible en écriture depuis le cluster
cd ~/work

## Créer un répertoire pour le TP : "Formation_Aln-SNP"
mkdir Formation_Aln-SNP ; cd Formation_Aln-SNP

## Créer les liens symboliques vers les données (fastq et référence) qui se trouvent dans 
## /save/formation/public_html/16_SGS-SNP/Data_Chr25-26/
ln -s /save/user/formation/public_html/16_SGS-SNP/Data_Chr25-26/G* .
ln -s /save/user/formation/public_html/16_SGS-SNP/Data_Chr25-26/S* .


################################################
## Exercice 1 : prise en main des données et du format FASTQ
################################################
## Prise en main des données 
## Se rendre sur le site du NCBI et en une requête rechercher les run SRR7062654 et SRR7062655
## Explorer les metadonnées : organisme ? échantillon ? séquenceur utilisé ? séquençage pairé ? ...

## Manipulation du format fastq
## Extraire et comparer les 5 premiers identifiants de chaque fastq du run SRR7062654
zgrep '^@SRR' SRR7062654_R1.fastq.gz | head -n 5
zgrep '^@SRR' SRR7062654_R2.fastq.gz | head -n 5

## Quel est le nombre de fragments dans le run SRR7062654 ?
zgrep -c '^@SRR' SRR7062654_R?.fastq.gz
##SRR7062654_R1.fastq.gz:993954
##SRR7062654_R2.fastq.gz:993954

## Combien de lectures contiennent un ou plusieurs « N » dans le run SRR7062654 ?
zgrep -c '^[ATCG]*N[ATCGN]*$' SRR7062654_R?.fastq.gz
##SRR7062654_R1.fastq.gz:2887
##SRR7062654_R2.fastq.gz:2778

## Analyse de la qualité du séquençage
## Lancer Fastqc sur les 2 runs et explorer les rapports générés
search_module fastqc
module load bioinfo/FastQC/0.12.1
for i in *fastq.gz ;  do echo "fastqc $i" >> 0_fastqc.jobs; done
mkdir LOGS
sarray -J 0_fastqc -e LOGS/%x_%j.err -o LOGS/%x_%j.out 0_fastqc.jobs


################################################
## Exercice 2 : alignement BWA
################################################
## Rechercher et charger BWA et samtools (commande module)
search_module bwa
module load bioinfo/bwa/0.7.17
search_module samtools
module load bioinfo/samtools/1.18

## Indexer la référence (bwa index)
bwa index Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa

## Lancer l'alignement en ajoutant les ReadGroup puis convertir en BAM et trier
## Une commande par couple de fastq à soumettre sur le cluster
\ls -1 *_R1.fastq.gz | perl -lne '$id=$_; $id=~s/\_.*//; $out=$_; $out=~s/\_R1.*//; $r2=$_; $r2=~s/\_R1\.f/\_R2\.f/; print "bwa mem -R \"\@RG\\tID:1\\tSM:$id\\tPL:illumina\\tLB:$id\\tPU:1\" -t4 Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa $_ $r2 | samtools sort - > $out.bam"' > 1_bwamem.jobs
sarray -J 1_bwamem -e LOGS/%x_%j.err -o LOGS/%x_%j.out -c 5 1_bwamem.jobs


################################################
## Excercice 3 : manipulation du format SAM/BAM - samtools
################################################
## Attention se connecter sur un noeud 
srun --pty bash

## Afficher l'aide de la commande samtools et parcourir les commandes disponibles
samtools --help 

## Le champ Flag, quels sont les différentes valeurs possible (et leur effectif) pour le run SRR7062654 ?
samtools view SRR7062654.bam | cut -f2 | sort -n |uniq -c

## Combien de lectures ont pour Cigar 90M ?
samtools view SRR7062654.bam | cut -f6 | grep -c '90M'
##1861309

## Combien de lectures sont parfaitement alignées ?
samtools view SRR7062654.bam | grep -c 'NM:i:0'
##1080285

## Comment expliquer cette différence ?
## M dans la Cigar = Match ou Mismatch

## Utiliser la commande calmd pour mettre en évidence les mismatches
samtools calmd -e SRR7062654.bam Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa | more

## Quelle est la valeur de qualité de mapping maximum et combien de lectures ont cette valeur ?
samtools view SRR7062654.bam | cut -f5 | sort -n |uniq -c
##1915573 60

## Convertir le fichier BAM du run SRR7062654 en SAM puis en BAM (une seule commande)
samtools view SRR7062654.bam | samtools view -b - > tmp_output.bam
## En cas d'erreur attention au header!
samtools view -h SRR7062654.bam | samtools view -b - > tmp_output.bam

## Trier le bam du run SRR7062654 par nom de reads et non par position (défaut)
samtools sort -n SRR7062654.bam -o SRR7062654_byName.bam

## Indexer les fichiers bam de votre répertoire de travail
\ls -1 *.bam | perl -lne 'print "samtools index $_"' > 2_index.jobs
sarray -J 2_index -e LOGS/%x_%j.err -o LOGS/%x_%j.out 2_index.jobs

## Utiliser la commande faidx de samtools pour extraire les 1 000 premières bases du chr25
samtools faidx Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa 25:1-1000

## Combien de lectures pour le run SRR7062654 s'alignent entre les coordonnées 1 et 1000 de la référence ?
samtools view SRR7062654.bam 25:1-1000 | wc -l
##74

## Calculer les statistiques d'alignement (flagstat, exécution sur le cluster)
\ls -1 *.bam | perl -lne 'print "samtools flagstat $_ > $_.flagstat"' > 3_flagstat.jobs
sarray -J 3_flagstat -e LOGS/%x_%j.err -o LOGS/%x_%j.out 3_flagstat.jobs

## En utilisant le champ « CIGAR », quelles sont les tailles de délétions et leurs effectifs ? 
samtools view SRR7062654.bam | perl -lane 'if($F[5]=~/(\d+)D/) {print $1}' | sort -n | uniq -c

## En utilisant le champ « MD », quel est le nombre de mismatch ?
samtools view SRR7062654.bam | perl -lne 'if($_=~/MD:Z:([0-9]+(([A-Z]|\^[A-Z]+)[0-9]+)*)/) {$md=$1; $md=~s/\^[ATGC]+//g; $md=~s/\d+//g; $t+=length($md);} END{print $t}' 
##1589125

################################################
# Excercice 4 : visualisation - IGV
################################################


################################################
# Excercice 5 : préprocessing et calling - GATK
################################################
## Marquager les duplicats (GATK MarkDuplicates) 
search_module gatk
module load bioinfo/GATK/4.4.0.0
module load devel/python/Python-3.11.1 devel/java/17.0.6
\ls *bam | perl -lne '$out=$_; $out=~s/\.bam//; print "gatk --java-options \"-Xms4000m -Xmx7g\" MarkDuplicates --INPUT $_ --METRICS_FILE $_.metrics --TMP_DIR . --ASSUME_SORT_ORDER coordinate --CREATE_INDEX true --OUTPUT $out.md.bam"' > 4_md.jobs
sarray -J 4_md -e LOGS/%x_%j.err -o LOGS/%x_%j.out --mem=10G 4_md.jobs

## Combien de lecture dupliquées dans le run SRR7062654
grep ^LIBRARY -A1 SRR7062654.bam.metrics | datamash transpose
##UNPAIRED_READ_DUPLICATES: 278
##READ_PAIR_DUPLICATES: 27580
samtools view -f 1024 SRR7062654.md.bam | wc -l
##55438 = 27580x2 + 278

## Recalibrer les bases (BQSR), étape 1 GATK BaseRecalibrator
\ls *md.bam | perl -lne '$out=$_; $out=~s/\.bam//; print "gatk --java-options -Xmx7g BaseRecalibrator -I $_ -O $out.recal.table --tmp-dir . -R Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa --known-sites Gallus_gallus_incl_consequences_chr25-26.vcf --verbosity INFO"' > 5_recal1.jobs
sarray -J 5_recal1 -e LOGS/%x_%j.err -o LOGS/%x_%j.out --mem=10G 5_recal1.jobs
## Erreur1 : remarque il faut l'index (.fai) de la référence 
samtools faidx Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa
## Erreur2 : remarque il faut le dictionaire (.dict) de la référence 
gatk --java-options "-Xms4000m -Xmx7g" CreateSequenceDictionary -R Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa

## Recalibrer les bases (BQSR), étape 2 GATK ApplyBQSR
\ls *md.bam | perl -lne '$out=$_; $out=~s/\.bam//; print "gatk --java-options -Xmx14g ApplyBQSR -R Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa --input $_ --output $out.recal.bam --bqsr-recal-file $out.recal.table"' > 6_recal2.jobs
sarray -J 6_recal2 -e LOGS/%x_%j.err -o LOGS/%x_%j.out --mem=15G 6_recal2.jobs

## A l'aide de samtools stats comparer avant / après recalibration
samtools stat SRR7062654.bam > SRR7062654.bam.stat
plot-bamstats SRR7062654.bam.stat -p SRR7062654.bam.stat
samtools stat SRR7062654.md.recal.bam > SRR7062654.md.recal.bam.stat
plot-bamstats SRR7062654.md.recal.bam.stat -p SRR7062654.md.recal.bam.stat

## Lancer le calling pour produire un G.VCF par échantillon (GATK HaplotypeCaller)
\ls *.md.recal.bam | perl -lne '$out=$_; $out=~s/\.bam//; print "gatk --java-options \"-Xmx10g -Xms1000m\" HaplotypeCaller -R Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa -I $_ -O $out.g.vcf -ERC GVCF"' > 7_HaplotypeCaller.jobs
sarray -J 7_haplotypecaller -e LOGS/%x_%j.err -o LOGS/%x_%j.out --mem=15G 7_HaplotypeCaller.jobs

## A partir des G.VCF individuels créer un G.VCF multisample 
\ls *.g.vcf | perl -ne 'BEGIN { print "gatk --java-options -Xmx10g CombineGVCFs -R Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa -O SRR.g.vcf"} chomp(); print " --variant $_"; END{print "\n"}' > 8_CombineGVCFs.jobs
sarray -J 8_CombineGVCFs -e LOGS/%x_%j.err -o LOGS/%x_%j.out --mem=15G 8_CombineGVCFs.jobs

## Faire le Calling GATK GenotypeGVCFs
cat > 9_GenotypeGVCFs.jobs
gatk --java-options "-Xmx10g" GenotypeGVCFs -R Gallus_gallus.Gallus_gallus-5.0.dna.toplevel_chr25-26.fa -V SRR.g.vcf -O SRR.vcf.gz
sarray -J 9_GenotypeGVCFs -e LOGS/%x_%j.err -o LOGS/%x_%j.out --mem=15G 9_GenotypeGVCFs.jobs


################################################
# Excercice 6 : format VCF annotation et filtre
################################################
## Repérer les différents champs du fichier VCF obtenu précédemment.
## Combien de variants contient le fichier VCF ?
zgrep -cv "^#"   SRR.vcf.gz
##82386

## Afficher l'aide de la commande SNPeff 
search_module snpeff
module load bioinfo/SnpEff/4.3T
snpEff -h

## Une annotation est-elle disponible pour notre génome d'intérêt ? Si oui, laquelle ?
snpEff databases | grep -i gallus
##Gallus_gallus-5.0.86

## Lancer l'annotation et repérer les ajouts dans le VCF (champ INFO)
snpEff ann Gallus_gallus-5.0.86 SRR.vcf.gz -s SRR-snpEff_summary.html > SRR-snpEff.vcf ; bgzip SRR-snpEff.vcf ; tabix -p vcf SRR-snpEff.vcf.gz

## Afficher l'aide de la commande SNPsift
java -Xmx4G -jar /usr/local/bioinfo/src/SnpEff/SnpEff-4.3T/snpEff/SnpSift.jar

## Ajouter l'information « snp connus » dans le champ ID du VCF
java -Xmx4G -jar /usr/local/bioinfo/src/SnpEff/SnpEff-4.3T/snpEff/SnpSift.jar annotate Gallus_gallus_incl_consequences_chr25-26.vcf SRR-snpEff.vcf.gz > SRR-snpEff-annotate.vcf ; bgzip SRR-snpEff-annotate.vcf ; tabix -p vcf SRR-snpEff-annotate.vcf.gz

## Combien de "nouveaux" variants dans le VCF
zgrep -v "^##" SRR-snpEff-annotate.vcf.gz | cut -f3 | grep -c '\.' 
##20146

## Filtre du vcf
## Combien de variants ont une qualité >= 100
java -Xmx4G -jar /usr/local/bioinfo/src/SnpEff/SnpEff-4.3T/snpEff/SnpSift.jar filter "QUAL>=100" -f SRR-snpEff-annotate.vcf.gz | grep -vc '^#'
##79213

## Combien de nouveaux SNP de qualité >= 200 sont homozygotes références pour SRR7062654 et homozygotes alternatifs pour SRR7062655 ?
java -Xmx4G -jar /usr/local/bioinfo/src/SnpEff/SnpEff-4.3T/snpEff/SnpSift.jar filter "ID == '' & isHom(GEN[SRR7062654]) & isHom(GEN[SRR7062655]) & isRef(GEN[SRR7062654]) & isVariant(GEN[SRR7062655]) & (REF=~'^[ATGC]$') & (ALT=~'^[ATGC]$') & (DP>20) & QUAL>=200" -f SRR-snpEff-annotate.vcf.gz | grep -cv '^#' 
##445

## Trouver la position du SNP nouveau ayant un effet "HIGH" et la qualité maximum
java -jar -Xmx4G /usr/local/bioinfo/src/SnpEff/SnpEff-4.3T/snpEff/SnpSift.jar filter "(ANN[*].IMPACT='HIGH') & (REF=~'^[ATGC]$') & (ALT=~'^[ATGC]$')" -f SRR-snpEff-annotate.vcf.gz | grep -v "^#"  |cut -f1-6 |sort -n -k6,6
## 25	1753972	.	T	G	965.5