Wiki

CR réunions

15/09/2022 : Natacha Nikolic - Sarah Maman
  • Objectif est de comparer les échantillons entre eux afin de mettre en évidence les régions plus ou moins conservées. Les échantillons étant anciens (Paléolithe jusqu'à l'âge moderne), les dégradations sont variables.
  • Si des zones partagées entre séquences sont mises en évidence, cela permettra de définir plus finement une couverture de séquençage pour une demande de budégitisation (WGS ou Capture).

19/10/2022: Demande d'extension de quota du /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/

23/11/2022: Point davancement. Présentation disponible en bas de page (avancement_novembre2022.odp).

1. Description du projet

1.1- Contacts

Natacha Nikolic <>
Stefanie Wagner <>
Aurelie Manicki <>
Patrick Jacques

Dans Sigenae: Sarah Maman, Mathieu Charles.

1.2- Contexte biologique du projet PaleoFish

Pas d'accès à MeeroDrop (closed), une partie des données transmises par Yannick MARIE, mail du 21/10/2022 "données brutes, fichiers qualité".
Séquençage faible profondeur sur des échantillons anciens, récupérés sur des sites archéologiques. Ces échantillons sont donc dégradés et peuvent contenir de l'ADN de l'hôte et de l'ADN exogène (contaminants). Ces échantillons ont été traités avec le pipeline PALEOMIX. Certains échantillons sont incomplets (problème de téléchargement).
Calcul du taux de recouvrement des échantillons sur plusieurs génomes de référence supposés afin de les classer. Ces traitements ont été réalisés sous Galaxy et ont été publiés.

L'analyse de la qualité des FASTQ (avec FATSQC report) et leur nettoyage (avec FastP) ont déjà été réalisés.

Pour la suite du projet, nous repartons des FASTQ disponibles : les FASTQ nettoyés avec FASTP.

Le calcul du taux de recouvrement entre les échantillons eux-même a pour objectif de mettre en évidence des régions plus conservées que d'autres afin de définir quelle profondeur de séquençage budgéter ?
De plus, ces calculs donneront des informations sur la dégradation des échantillons.

Noms des échantillons:
_com = génome total (nucléaire et mitochondrial).
_mt = uniquement le génome mitochondrial.

TODO Fastq Screen ??

1.3- Question

Les données brutes (FastQ) ont été nettoyées sous Galaxy et Paleomix. Les alignements ont été fait sur Anguilla anguilla, Anguilla rostrata, Salmo salar, et Salmo trutta.
Nous avons besoin des % de couverture et % de génome partagé entre nos échantillons et entre les échantillons et génome de référence.

Informations sur les échantillons :
Les résultats de nos échantillons passés via Paleomix sont hébergés chez vous « Project/ECOBIOP_ADNa/Data_PaleoFish (comb = alignement sur génome complet donc nucléiare et mitochondrial) :
Parmi ces résultats nous avons les fichiers Coverage, dephs, BAM etc. :
Nous savons que 2 échantillons (L5-1 et L9-5) ne sont probablement pas du saumon ou de l’anguille, puisque avec Galaxy en faisant un BWA nous avions vu que le premier s’alignait avec Hucho hucho et le second avec l’homme.
Pour le moment, ce n’est pas la peine de regarder ces 2 échantillons.

Les informations récupérées avec Galaxy (après nettoyage par fastp, info qualité par fastQC et alignement par BWA) sont également sur notre espace de travail. Le % de couverture semblait variable mais nous n’avions pas l’info sur le % de génome partagé entre échantillon.

2. Pipeline pour Anguilla anguilla

Répertoire de travail

/work/project/sigenae/sarah/PALEOFISH/Reference_Genome_Whole_genome/Poissons

Indexation des génomes de référence

Anguilla_anguilla

Vérification des FASTQ

Question 1 :
Certains FASTQ Anguilla anguilla sont corrompus et d'autres sont manquants:

Pour les foward:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae $ for filename in /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/*.fastq.gz;do
>  gunzip -t $filename/forward.fastqsanger.gz ; [ $? == 0 ] || echo pb $filename/forward.fastqsanger.gz
> done

gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_10_R9_54_e_S122.fastq.gz/forward.fastqsanger.gz: unexpected end of file
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_10_R9_54_e_S122.fastq.gz/forward.fastqsanger.gz

gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_14_SCO_Aa_227_S126.fastq.gz/forward.fastqsanger.gz: unexpected end of file
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_14_SCO_Aa_227_S126.fastq.gz/forward.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_15_SCO_Aa_229_S127.fastq.gz/forward.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_15_SCO_Aa_229_S127.fastq.gz/forward.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_10_SCO_2112_S138.fastq.gz/forward.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_10_SCO_2112_S138.fastq.gz/forward.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_4_SCO_1451_S132.fastq.gz/forward.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_4_SCO_1451_S132.fastq.gz/forward.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_5_SCO_2_S133.fastq.gz/forward.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_5_SCO_2_S133.fastq.gz/forward.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_8_SCO_519_S136.fastq.gz/forward.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_8_SCO_519_S136.fastq.gz/forward.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_9_SCO_1492_S137.fastq.gz/forward.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_9_SCO_1492_S137.fastq.gz/forward.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_11_Aa_Sou266_248_1_S91.fastq.gz/forward.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_11_Aa_Sou266_248_1_S91.fastq.gz/forward.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_12_Aa_Sou266_794_1_S92.fastq.gz/forward.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_12_Aa_Sou266_794_1_S92.fastq.gz/forward.fastqsanger.gz

gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_2_Aa_Sou266_943_1_S82.fastq.gz/forward.fastqsanger.gz: unexpected end of file
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_2_Aa_Sou266_943_1_S82.fastq.gz/forward.fastqsanger.gz

Pour les reverse:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae $ for filename in /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/*.fastq.gz;do
> gunzip -t $filename/reverse.fastqsanger.gz ; [ $? == 0 ] || echo pb $filename/reverse.fastqsanger.gz
> done

gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_10_R9_54_e_S122.fastq.gz/reverse.fastqsanger.gz: unexpected end of file
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_10_R9_54_e_S122.fastq.gz/reverse.fastqsanger.gz

gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_14_SCO_Aa_227_S126.fastq.gz/reverse.fastqsanger.gz: unexpected end of file
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_14_SCO_Aa_227_S126.fastq.gz/reverse.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_15_SCO_Aa_229_S127.fastq.gz/reverse.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_15_SCO_Aa_229_S127.fastq.gz/reverse.fastqsanger.gz

gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_8_R9_54_a_S120.fastq.gz/reverse.fastqsanger.gz: unexpected end of file
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L10_8_R9_54_a_S120.fastq.gz/reverse.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_10_SCO_2112_S138.fastq.gz/reverse.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_10_SCO_2112_S138.fastq.gz/reverse.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_4_SCO_1451_S132.fastq.gz/reverse.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_4_SCO_1451_S132.fastq.gz/reverse.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_5_SCO_2_S133.fastq.gz/reverse.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_5_SCO_2_S133.fastq.gz/reverse.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_8_SCO_519_S136.fastq.gz/reverse.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_8_SCO_519_S136.fastq.gz/reverse.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_9_SCO_1492_S137.fastq.gz/reverse.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L11_9_SCO_1492_S137.fastq.gz/reverse.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_11_Aa_Sou266_248_1_S91.fastq.gz/reverse.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_11_Aa_Sou266_248_1_S91.fastq.gz/reverse.fastqsanger.gz
gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_12_Aa_Sou266_794_1_S92.fastq.gz/reverse.fastqsanger.gz: No such file or directory
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_12_Aa_Sou266_794_1_S92.fastq.gz/reverse.fastqsanger.gz

gzip: /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_2_Aa_Sou266_943_1_S82.fastq.gz/reverse.fastqsanger.gz: unexpected end of file
pb /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/L8_2_Aa_Sou266_943_1_S82.fastq.gz/reverse.fastqsanger.gz

Question 2 :
Les 69 FASTQ Salmonidae dans /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Salmonidae/, mais lequels sont Salmo trutta et Salmo salar ? J'aurai besoin de les distinguer pour travailler sur la bonne référence.
Pas de fichiers manquants, pas de fichiers foward ni reverse corrompu.

Après mis à jour des FASTQ par Natacha, atres points à vérifier: bien le même nombre de r1 et r2, et donnés dans le même ordre ?

sigenae@genologin1 /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/EG10.fastq $ for filename in /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/*.fastq.gz;do
> grep @ reverse.fastqsanger | wc -l
> done
740333
740333
740333
...
Idem pour forward.
  • Nettoyage du répertoire AUD10199-1.fastq pour ne garder que AUD10199-1.fastq/forward.fastqsanger AUD10199-1.fastq/reverse.fastqsanger:
    sigenae@genologin1 /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae $ ls AUD10199-1.fastq/*
    AUD10199-1.fastq/forward.fastqsanger  AUD10199-1.fastq/reverse.fastqsanger
    
    AUD10199-1.fastq/AUD10199-1.fastq:
    
    AUD10199-1.fastq/AUD10310-2.fastq:
    
    AUD10199-1.fastq/AUD10367-1.fastq:
    
    AUD10199-1.fastq/AUD10487-1.fastq:
    forward.fastqsanger  reverse.fastqsanger
    
    AUD10199-1.fastq/AUD10671-2.fastq:
    
    AUD10199-1.fastq/AUD10857-2.fastq:
    forward.fastqsanger
    
    AUD10199-1.fastq/AUD10996-1.fastq:
    forward.fastqsanger  reverse.fastqsanger
    
    AUD10199-1.fastq/AUD11096-1.fastq:
    
    AUD10199-1.fastq/AUD11719-1.fastq:
    
    AUD10199-1.fastq/AUD11719-2.fastq:
    
    AUD10199-1.fastq/AUD11747-9.fastq:
    
    AUD10199-1.fastq/AUD11782-2.fastq:
    
    AUD10199-1.fastq/AUD20457-2.fastq:
    
    AUD10199-1.fastq/CD2:
    fastp on collection 5_ Paired-end output/
    
    AUD10199-1.fastq/EG10.fastq:
    
  • Le répertoire CD2 était présent 2 fois:

ls /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/CD2/fastp\ on\ collection\ 5_\ Paired-end\ output/CD2.fastq.gz/
forward.fastqsanger.gz reverse.fastqsanger.gz
sigenae@genologin1 /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp $ ls /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/CD2.fastq.gz/
--> j'ai récupèré les fichiers forward et reverse de
CD2/fastp\ on\ collection\ 5_\ Paired-end\ output/CD2.fastq.gz/
pour les mettre dans:
CD2.fastq.gz

Objectif

Distinguons:
  • Taux de couverture : profondeur de séquençage = depth = nombre de lectures.
  • Taux de recouvrement : overlapping = breadth = chevauchement des reads pour s'aligner sur la référence, pourcentage en longeur du génome qui se trouve aligné avec des reads, nombre de bases sur la référence.

Taux de recouvrement sur la référence:
bedtools coverage -a reference.VCF (description des chromosomes) -b sample1.bed sample2.bed ... -hist

Taux de recouvrement d'un échantillon par rapport à un autre (s1 sur s2): voir les étapes ci-dessous

Indexation de la référence

bwa index Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta

(base) smaman@node140 /work/project/sigenae/sarah/PALEOFISH/Reference_Genome_Whole_genome/Poissons $ bwa index  Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta
[bwa_index] Pack FASTA... 9.22 sec
[bwa_index] Construct BWT for the packed sequence...
[BWTIncCreate] textLength=1958058486, availableWord=149775832
[BWTIncConstructFromPacked] 10 iterations done. 99999990 characters processed.
...
[BWTIncConstructFromPacked] 260 iterations done. 1950183782 characters processed.
[bwt_gen] Finished constructing BWT in 264 iterations.
[bwa_index] 1152.95 seconds elapse.
[bwa_index] Update BWT... 6.12 sec
[bwa_index] Pack forward-only FASTA... 4.96 sec
[bwa_index] Construct SA from BWT and Occ... 408.09 sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa index Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta
[main] Real time: 1582.803 sec; CPU: 1581.352 sec

Référence indexée:

(base) smaman@node140 /work/project/sigenae/sarah/PALEOFISH/Reference_Genome_Whole_genome/Poissons $ ls -ltrah
-rw-r--r-- 1 smaman  SIGENAE 934M 16 sept. 16:56 Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta.bwt
-rw-r--r-- 1 smaman  SIGENAE 234M 16 sept. 16:56 Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta.pac
-rw-r--r-- 1 smaman  SIGENAE 6,4K 16 sept. 16:56 Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta.ann
-rw-r--r-- 1 smaman  SIGENAE  10K 16 sept. 16:56 Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta.amb
drwxrwsr-x 3 sigenae SIGENAE 8,0K 16 sept. 17:03 .
-rw-r--r-- 1 smaman  SIGENAE 467M 16 sept. 17:03 Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta.sa

Aligner sur la référence

s1.fastq.gz -> s1.bam

bwa aln sampe est plus précis mais mappe moins. En effet, les positions alternatives ne sont pas mappées. S'il y a trop de distances entre les reads et le génome de référence, les reads ne sont pas mappées, il n'y a pas d'alignement local.
Il est préférable d'utiliser bwa mem pour plus de mapping. S. Wagner a mentionné le fait qu'avec l'ADN ancien il est préférable d'utiliser aln que mem à cause des nucléotides manquants comblés par mem. Sarah pense que ce n'est pas un problème car on supprime après les alignements secondaires mais à vérifier.
L'alignement local permet les alignements secondaires donc forcera le recouvrement. Ensuite, il sera donc nécessaire de filtrer les reads non mappées et les alignements secondaires:

Préparation du répertoire RESULTATS/:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ mkdir /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS/

Distinguer les échantillons mâle/femelle car la référence n'est pas la même (liste demandée par mail à Natacha le 7/11):

Préparation du fichier de lancement 1-mapp.sh (#bwa mem "index_prefix" "s1_R1.fq" "s1_R2.fq" > s1.sam):

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ for filename in /work/project/sigenae/sarah/PALEOFISH/Clean_data_after_fastp/Anguillidae/*.fastq.gz; do echo " module purge && module load bioinfo/bwa-0.7.17  &&  bwa mem $ref $filename/forward.fastqsanger.gz  $filename/reverse.fastqsanger.gz > /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS/$(basename $filename).$i.sam &&  cd /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS/ && ln -s $(basename $filename).$i.sam  $i.sam "; i=$(( $i + 1 )); done

Visualisation du fichier de lancement:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ more 1-mapp.sh 
 ....

Lancement du job:

cd  /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/
sarray -c 4 --mem=10G 1-mapp.sh
Submitted batch job 

et

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ sarray -c 4 --mem=10G  1-mapp-suite.sh
Submitted batch job 37376508

Filtre

samtools view -F 256 -F 4 sam
-F : exclude
256: flag des secondary alignements
4: flag des reads non mappées
https://broadinstitute.github.io/picard/explain-flags.html

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ more 2-filter-convert.sh 
#!/bin/bash
#SBATCH -J paleofish
#SBATCH -p unlimitq 
#SBATCH --mem=20G

module purge
module load bioinfo/samtools-1.16.1
module load bioinfo/bedtools-2.27.1

# Filter + Convertir le BAM en BED: 
cd  /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS/
for ((i=1; i<=60; i++)); do
    samtools view -Sb -F 256 -F 4 $i.sam | samtools sort -@ 4 -m 1G  - | bamToBed -i - > $i.bed
done

Lancement du job:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ sbatch 2-filter-convert.sh
Submitted batch job 37379842  (de 1 à 48)
et
(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ sbatch 2-filter-convert-suite.sh
Submitted batch job 37395053 (de 49 à 60)

Regrouper les intervals avec bedtools merge

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ more 3-merge-intervals.sh 
#!/bin/bash
#SBATCH -J paleofish
#SBATCH -p unlimitq 
#SBATCH --mem=20G

module purge
module load bioinfo/bedtools-2.27.1

# bedtools merge -i s1.bed > s1.merge.bed 
cd  /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS/
for ((i=1; i<=60; i++)); do
    bedtools merge -i $i.bed > $i.merge.bed 
done

Lancement du job:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ sbatch  3-merge-intervals.sh 
Submitted batch job 37402064

Récupérer le nombre de bases communes entre s1 et s2 avec bedtools coverage

#https://bedtools.readthedocs.io/en/latest/content/tools/coverage.html
bedtools coverage -a s1.merge.bed -b s2.merge.bed

Voici l'exemple donné par les bedtools:


$ cat A.bed
chr1  0   100
chr1  100 200
chr2  0   100

$ cat B.bed
chr1  10  20
chr1  20  30
chr1  30  40
chr1  100 200

$ bedtools coverage -a A.bed -b B.bed
chr1  0   100  3  30  100 0.3000000
chr1  100 200  1  100 100 1.0000000
chr2  0   100  0  0   100 0.0000000

Colonne 1: Nom du chromosome
Colonne 2: 0 : position start dans A.bed (chr1  0   100)
Colonne 3: 100 : position stop dans A.bed (chr1  0   100)
Colonne 4: le profondeur = la position 0-100 chr1 de A est vu 3 fois dans B : chr1 10-20, chr1 20-30 et chr1 30-40, à ne pas prendre en compte car nos échantillons sont mergés.
Colonne 5: recouvrement = 30 nucleotides sont communs de 10 à 40 (40-10=30)
Colonne 6: longeur dans A : 100 bases 
Colonne 7: taux de recouvrement = 30/100 = col 5/ col 6 = 0.3

Traitement des données Anguilla anguilla:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ more 4-bedtools-coverage.sh 
#!/bin/bash
#SBATCH -J paleofish
#SBATCH -p unlimitq 
#SBATCH --mem=15G

module purge
module load bioinfo/bedtools-2.27.1

cd /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS/;
for ((i=1; i<=59; i++)); do
    for ((j=i; j<=60; j++)); do
        bedtools coverage -a $i.merge.bed -b $j.merge.bed > coverage_$i-$j.out
    done    
done

Lancement du job:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ sbatch  4-bedtools-coverage.sh 
Submitted batch job 37511889

Calcul du taux de recouvrement

taux de recouvremement = nb tot bases (5ieme colonne)/ (long chromosomes +scaffolds)
avec

1/ Long chromosomes +scaffolds

#https://www.bioinformatics.nl/cgi-bin/emboss/help/infoseq
(infoseq -sequence /bank/ebi/ensembl/ensembl_sus_scrofa_genome/current/fasta/ensembl_sus_scrofa_genome -length -only -noheading 2>/dev/null | perl -lane '$n+=$F0;END{print $n;}'
ou)

srun --mem 10G --pty bash 
module load system/datamash-1.3
module load bioinfo/EMBOSS-6.6.0
infoseq -sequence /work/project/sigenae/sarah/PALEOFISH/Reference_Genome_Whole_genome/Poissons/Anguilla_anguilla_GCA_013347855.1_fAngAng1.pri_genomic.fasta -length -only -noheading 2>/dev/null| sed -e 's/\s//g'|datamash sum 1 

Résultat:
Long chromosomes +scaffolds = 979029243

2/ Nb tot bases (5ieme colonne)

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH $ search_module datamash
system/datamash-1.3
(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH $ module load system/datamash-1.3
(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH $ datamash sum 5  coverage_$i-$j.out

Préparation du job:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ more 5-calculate-coverage-by-sample.sh 
#!/bin/bash
#SBATCH -J paleofish
#SBATCH -p unlimitq 
#SBATCH --mem=15G

module purge
module load system/datamash-1.3

cd /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS/;
for ((i=1; i<=59; i++)); do
    for ((j=i; j<=60; j++)); do
        datamash -s sum 5  < coverage_$i-$j.out > coverage_$i-$j.res
    done    
done

Lancement du job:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ sbatch 5-calculate-coverage-by-sample.sh 
Submitted batch job 37529887

3/ Mise en forme de la matrice

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ more MATRICE.sh
#!/bin/bash

nb=$(cat $file)
for ((i=1; i<=59; i++)); do
    for ((j=i; j<=60; j++)); do
        file="coverage_$i-$j.res" 
        echo "$i; $j; $(cat $file)" >> MATRICE.txt
    done    
done

Résultat:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ ./MATRICE.sh 

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ sbatch  MATRICE.sh
Submitted batch job 37545957

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ more MATRICE.txt 
1; 1; 2909603
1; 2; 105370
1; 3; 106808
1; 4; 125690
1; 5; 62378
1; 6; 112123
1; 7; 59329
1; 8; 74597
1; 9; 291343
1; 10; 83517
1; 11; 148013
1; 12; 356791

...

Ajout du calcul des ratios:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ more calculate_for_matrice.sh 
#!/bin/bash

# ./calculate_for_matrice.sh  MATRICE_head100.txt
#  ./calculate_for_matrice.sh MATRICE.txt

for line in $(cat $1)
do
    IFS='$' read -ra ADDR <<< "$line" 
    #echo "line= $line" 
    for i in {1..5}; 
    do  
        if [[ "$line" =~ ^"$i","$i",.* ]]; then
            IFS=',' read -ra nb <<< "$line" 
            den="${nb[2]}" 
        fi
        if [[ "$line" =~ ^"$i",.* ]]; then
            IFS=',' read -ra n <<< "$line" 
            num="${n[2]}" 
            ratio=$(echo "scale=4 ; $num / $den" | bc)
            newline="$line,$ratio" 
            echo $newline >> MATRICE_withRatio.csv
        fi
   done
done

Visualisation:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ head -n 20  MATRICE_withRatio.csv 
1,1,2909603,1.0000
1,2,105370,.0362
1,3,106808,.0367
1,4,125690,.0431
1,5,62378,.0214
1,6,112123,.0385
...

Proposition de représentation des résultats avec la librairie seaborn: ratio nb nucléotides de l'échantillon / nb nucléotides de l'échantillon de référence (voir le fichier MATRICE_RATIOS.csv). Le tableau de correspondance est disponible en bas de page : tableau_correspondance_paleofish_anguilla_anguilla.ods

paleofish_ratios.png

Pour générer la partie gauche de la heatmap, relancer les dernières étapes mais dans /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/:

(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ cp 4-bedtools-coverage.sh  4-bedtools-coverage_BIS.sh 
(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ cp  5-calculate-coverage-by-sample.sh  5-calculate-coverage-by-sample_BIS.sh 
(base) smaman@genologin1 /work/project/sigenae/sarah/PALEOFISH/recouvrement_entre_samples_Anguilla_anguilla $ cp RESULTATS/MATRICE.sh RESULTATS/MATRICE_BIS.sh

Modification des boucles:

for ((i=1; i<=59; i++)); do
    for ((j=i; j<=60; j++)); do

en:
for ((i=1; i<=60; i++)); do
    for ((j=1; j<=60; j++)); do

Désormais, les jobs sont lancés sur /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla $ sbatch 4-bedtools-coverage_BIS.sh
Submitted batch job 40520707
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla $ sbatch -d afterok:40520707 5-calculate-coverage-by-sample_BIS.sh
Submitted batch job 40520708
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla $ cd RESULTATS/
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ sbatch -d afterok:40520708 MATRICE_BIS.sh
Submitted batch job 40520710

Puis utiliser MATRICE_BIS.txt pour lancer calculate_for_matrice.sh

/work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ more MATRICE_BIS.txt 
1; 1; 2909603
2; 1; 105370
2; 2; 12770347
$ sed -i 's/;/,/g' MATRICE_BIS.txt

smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ python3 MATRICE_BIS.txt

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS $ head MATRICE_ratio.csv
1,1,2909603,1.0000
1,2,105370,.0362
1,3,106808,.0367
1,4,125690,.0431
1,5,62378,.0214
1,6,112123,.0385
1,7,59329,.0203
1,8,74597,.0256
1,9,291343,.1001
1,10,83517,.0287

Code pour générer cette seconde heatmap

Anguilla_anguilla_gauche.png

3. Pipeline pour le saumon Salmonidae Salmo salar

Répertoire de travail

/work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/

salmo salar:

/work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/FASTQ/

Indexation des génomes de référence Salmo salar

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/Ref $ ln -s  /work/project/sigenae/sarah/PALEOFISH/Reference_Genome_Whole_genome/Poissons/Salmo_Salar_Male_GCF_905237065.1_Ssal_v3.1_genomic.fna.gz  .
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/Ref $ ln -s /work/project/sigenae/sarah/PALEOFISH/Reference_Genome_Whole_genome/Poissons/Salmo_salar_Female_AtlanticSalmon_WholeGenome_GCF_000233375.1_ICSASG_v2_genomic.fna.gz].fasta.gz  .

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/REF $ more index_Salmo_salar_Female.sh
#!/bin/bash
#SBATCH -p unlimitq
#SBATCH --mem=100G
module load bioinfo/bwa-0.7.17
cd /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/REF/
bwa index Salmo_salar_Female_AtlanticSalmon_WholeGenome_GCF_000233375.1_ICSASG_v2_genomic.fna.gz].fasta.gz

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/REF $ more index_Salmo_salar_Male.sh
#!/bin/bash
#SBATCH -p unlimitq
#SBATCH --mem=100G
module load bioinfo/bwa-0.7.17
cd /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/REF/
bwa index Salmo_Salar_Male_GCF_905237065.1_Ssal_v3.1_genomic.fna.gz

Index fasta female: 40779585
Index fasta male: 40809373

Vérification des FASTQ

Pour les forward:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/FASTQ $ for filename in *;do
> gunzip -t $filename/forward.fastqsanger.gz ; [ $? == 0 ] || echo pb $filename/forward.fastqsanger.gz
> done
gzip: BC4.fastq/forward.fastqsanger.gz: No such file or directory
pb BC4.fastq/forward.fastqsanger.gz
gzip: MPHB2/forward.fastqsanger.gz: No such file or directory
pb MPHB2/forward.fastqsanger.gz
gzip: MPHB3/forward.fastqsanger.gz: No such file or directory
pb MPHB3/forward.fastqsanger.gz
gzip: MPHB4/forward.fastqsanger.gz: No such file or directory
pb MPHB4/forward.fastqsanger.gz
gzip: OLG1.fastq/forward.fastqsanger.gz: No such file or directory
pb OLG1.fastq/forward.fastqsanger.gz

Pour les reverse:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/FASTQ $ for filename in *;do gunzip -t $filename/forward.fastqsanger.gz ; [ $? == 0 ] || echo pb $filename/reverse.fastqsanger.gz; done
gzip: BC4.fastq/forward.fastqsanger.gz: No such file or directory
pb BC4.fastq/reverse.fastqsanger.gz
gzip: MPHB2/forward.fastqsanger.gz: No such file or directory
pb MPHB2/reverse.fastqsanger.gz
gzip: MPHB3/forward.fastqsanger.gz: No such file or directory
pb MPHB3/reverse.fastqsanger.gz
gzip: MPHB4/forward.fastqsanger.gz: No such file or directory
pb MPHB4/reverse.fastqsanger.gz
gzip: OLG1.fastq/forward.fastqsanger.gz: No such file or directory
pb OLG1.fastq/reverse.fastqsanger.gz

Compresser les fastq non compressés, exemple:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/FASTQ/MPHB2 $ gzip forward.fastqsanger

Pas de fichiers manquants, pas de fichiers foward ni reverse corrompu.

Vérifier le même nombre de r1 et r2, et donnés dans le même ordre ?

sample forward reverse
BC4 3913038 3913038
L4_11_MD2_S27 1133821 1133821
L4_8_MB4_S24 1040332 1040332
L4_9_MD12_S25 923045 923045
L5_12_MD9_S44 1260731 1260731
L5_13_BC11_S45 1238316 1238316
L5_14_BC12_S46 1418314 1418314
L5_15_BC13_S47 1023126 1023126
L5_2_MD7_S34 1283879 1283879
L5_3_MD8_S35 1262435 1262435
L5_5_MD10_S37 1606892 1606892
L5_8_MD16_S40 1180661 1180661
L9_1_26_52_S97 1080871 1080871
L9_2_2777_17_S98 884229 884229
MPHB2 5958286 5958286
MPHB3 2866653 2866653
MPHB4 3393272 3393272
OLG1 934034 934034
$ zgrep @   reverse.fastqsanger.gz  | wc -l

Aligner sur la référence, filtrer et regrouper

s1.fastq.gz -> s1.bam

Préparation du répertoire RESULTATS/:
Préparation du fichier de lancement 1-mapp.sh (#bwa mem "index_prefix" "s1_R1.fq" "s1_R2.fq" > s1.sam):

Attention, comme nous avons deux références (male et femelle) et pas de lien entre le nom de l'échantillon et le sexe, le mapping de chaque échantillon est donc lancé une fois sur la référence mâle et une fois sur la référence femelle.
Pour chaque échantillon, nous aurons un female.sam et un male.sam.

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/ $ for filename in /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/FASTQ/*.fastq.gz; do echo " module purge && module load bioinfo/bwa-0.7.17  &&  bwa mem $ref $filename/forward.fastqsanger.gz  $filename/reverse.fastqsanger.gz > /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/RESULTATS/$(basename $filename).$i.sam &&  cd /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/RESULTATS/ && ln -s $(basename $filename).$i.sam  $i.sam "; i=$(( $i + 1 )); done

Lancement des jobs:

cd  /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar $ sarray -c 4 --mem=10G 1-mapp.sh
Submitted batch job 40810250 / COMPLETED
Le numéro 14 était manquant:
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar $  sarray -c 4 --mem=10G 1-mapp_num14.sh
Submitted batch job 40813416 / OK
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar $ sbatch  2-filter-convert.sh 
Submitted batch job 40844311 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar $ sbatch -d afterok:40844311 3-merge-intervals.sh 
Submitted batch job 40844312 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar $ sbatch -d afterok:40844312 4-bedtools-coverage.sh 
Submitted batch job 40844313 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar $ sbatch -d afterok:40844313 5-calculate-coverage-by-sample.sh 
Submitted batch job 40844314 / COMPLETED

Voici la correspondance entre le nom de l’échantillon et son numéro:

Male= numéro impaire, Femelle= numéro paire.
Nom de l'échantillon numero (male/femelle)
L4_11_MD2_S27 1.male =1 1.femelle=2
L4_8_MB4_S24 2.male =3 2.femelle=4
L4_9_MD12_S25 3.male =5 3.femelle=6
L5_12_MD9_S44 4.male =7 4.femelle=8
L5_13_BC11_S45 5.male =9 5.femelle=10
L5_14_BC12_S46 6.male =11 6.femelle=12
L5_15_BC13_S47 7.male =13 7.femelle=14
L5_2_MD7_S34 8.male =15 8.femelle=16
L5_3_MD8_S35 9.male =17 9.femelle=18
L5_5_MD10_S37 10.male =19 10.femelle=20
L5_8_MD16_S40 11.male =21 11.femelle=22
L9_1_26_52_S97 12.male =23 12.femelle=24
L9_2_2777_17_S98 13.male =25 13.femelle=26
MPHB2 14.male =27 14.femelle=28
MPHB3 15.male =29 15.femelle=30
MPHB4 16.male =31 16.femelle=32
OLG1 17.male =33 17.femelle=34
BC4 18.male =35 18.femelle=36

Salmo_salar.png

En retirant les rations les plus élevés (dont la diagonale), alors le damier est plus visible. Pour chaque échantillon, il semblerait que la qualité du génome (mâle ou femelle ?) soit mieux conservée. Il s'agit plus d'une conservation de génome que d'une distinction possible mâle/femelle.

Salmo_salar_damier.png

4. Pipeline pour le saumon Salmonidae Salmo trutta

Répertoire de travail

/work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/

Récupération et vérification des FASTQ:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/FASTQ $ cp -r /work/project/ECOBIOP_ADNa/Data_PaleoFish/Clean_data_after_fastp/Salmonidae/Salmo\ trutta/* .

Vérification des fichiers forward, reverse et du nombre de séquences pour Salmo trutta

Indexation de la référence

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/Ref $ cp /work/project/ECOBIOP_ADNa/Data_PaleoFish/Reference_Genome_Whole_genome/Poissons/SalmoTrutta_WholeGenome_GCF_901001165.1_fSalTru1.1_genomic.fna.gz].fasta.gz .
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/Ref $ mv SalmoTrutta_WholeGenome_GCF_901001165.1_fSalTru1.1_genomic.fna.gz].fasta.gz  Salmo_trutta.fasta.gz

Lancement de l'indexation pour Salmo trutta

Mapping, filtre, regroupement et matrice

Préparation du fichier 1_mapp.sh pour Salmo trutta

Lancement du mapping:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta $ sarray -c 4 --mem=10G 1-mapp.sh
Submitted batch job 40808996 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta $ sbatch 2-filter-convert.sh
Submitted batch job 40811308 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta $ sbatch -d afterok:40811308 3-merge-intervals.sh
Submitted batch job 40811310 / COMPLETED

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta $ sbatch 4-bedtools-coverage.sh 
Submitted batch job 40814148 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta $ sbatch -d afterok:40814148 5-calculate-coverage-by-sample.sh 
Submitted batch job 40814151 / COMPLETED

# Avec l'ensemble des comparaisons:
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta $ sbatch  4-bedtools-coverage.sh 
Submitted batch job 40844326 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta $ sbatch -d afterok:40844326 5-calculate-coverage-by-sample.sh 
Submitted batch job 40844328 / COMPLETED

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/RESULTATS $ sbatch MATRICE.sh
Submitted batch job 40844315 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/RESULTATS $ sed -i 's/; /,/g' MATRICE.txt
# Ajouter un header:
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/RESULTATS $ head MATRICE_withHeader.txt
sample1,sample2,nbnt
1,1,224865290
1,2,33238485

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/RESULTATS $ python3 calculate_ratio.py 
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/RESULTATS $ more MATRICE_ratios.csv
sample1,sample2,ratio
1 , 1 , 1.0
1 , 2 , 0.147815098541887
1 , 3 , 0.006690472326787295
1 , 4 , 0.010923379948946323
1 , 5 , 0.040903418219859544
1 , 6 , 0.03362676382824579
1 , 7 , 0.03248339483608164
1 , 8 , 0.0388004480371337
1 , 9 , 0.02139206099794237

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/RESULTATS $ tail -n 20  MATRICE_ratios.csv
34 , 31 , 0.028311983376624984
34 , 32 , 0.00409985223313063
34 , 33 , 0.2134600777232833
34 , 34 , 1.0

Correspondance entre le nom et le numéro de l'échantillon:
Nom Numéro
BC1 1
BC2 2
EG13 3
EG4 4
L10_1_BC5_S113 5
L10_2_BC6_S114 6
L10_3_BC7_S115 7
L10_4_BC8_S116 8
L10_5_BC9_S117 9
L11_1_SCO_98_S129 10
L11_3_SCO_114_S131 11
L3_10_MB13_S10 12
L3_12_MB6_S12 13
L3_13_MB11_S13 14
L3_14_MB10_S14 15
L3_1_MB15_S1 16
L3_5_SM3_S5 17
L3_6_MB3_S6 18
L4_10_BC14_S26 19
L4_12_MD1_S28 20
L4_15_BC10_S31 21
L4_3_MB16_S19 22
L4_6_MB2_S22 23
L4_7_MB1_S23 24
L5_9_MD6_S41 25
L6_11_SH2000_99_394_g_S59 26
L6_13_SH2000_99_394_f_S61 27
L6_5_SH1979_4_6843_d_S53 28
L6_8_SH2000_99_394_e_S56 29
L9_15_SH2000_99_390_c_S111 30
L9_7_SH2000_99_394_c_S103 31
SH2 32
TP1 33
TP2 34

Générer la matrice Salmo trutta

Salmo_trutta.png

5. Pipeline pour le saumon Salmonidae Salmo sp

Répertoire de travail

Récupération des FASTQ:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/FASTQ $ ls /work/project/ECOBIOP_ADNa/Data_PaleoFish/Clean_data_after_fastp/Salmonidae/Salmo\ sp/Salmo_sp_impossible_daller_plusLoin/
L10_11_SCO_1338_S123.fastq.gz  L11_7_SCO_1354_S135.fastq.gz  L4_13_MB5_S29.fastq.gz  L4_4_MB9_S20.fastq.gz    L5_4_MD4_S36.fastq.gz              L6_4_SH1979_4_6843_c_S52.fastq.gz      L9_14_SH1979_4_6843_a_S110.fastq.gz
L10_12_SCO_2213_S124.fastq.gz  L3_11_MB18_S11.fastq.gz       L4_14_SM2_S30.fastq.gz  L4_5_MD5_S21.fastq.gz    L5_6_MB8_S38.fastq.gz              L6_6_SH2001_106-11_b_S54.fastq.gz      L9_4_SM1_S100.fastq.gz
L11_2_SCO_81_S130.fastq.gz     L3_15_MB12_S15.fastq.gz       L4_1_MB14_S17.fastq.gz  L5_10_MD11_S42.fastq.gz  L5_7_MD3_S39.fastq.gz              L8_15_HTMK99_XXII_0_4287_S95.fastq.gz  L9_6_SH1979_4_6843_b_S102.fastq.gz
L11_6_SCO_116_S134.fastq.gz    L3_9_MB17_S9.fastq.gz         L4_2_MD13_S18.fastq.gz  L5_11_MD14_S43.fastq.gz  L6_3_SH1979_4_6843_e_S51.fastq.gz  L9_13_SH2000_99_390_b_S109.fastq.gz
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/FASTQ $ cp -r /work/project/ECOBIOP_ADNa/Data_PaleoFish/Clean_data_after_fastp/Salmonidae/Salmo\ sp/Salmo_sp_impossible_daller_plusLoin/* .

Vérification des FASTQ: L'ensemble des paires de fichier (forward, reverse) sont présents et compressés.

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/FASTQ $ for filename in *;do
> gunzip -t $filename/forward.fastqsanger.gz ; [ $? == 0 ] || echo pb $filename/forward.fastqsanger.gz
> done
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/FASTQ $ for filename in *;do gunzip -t $filename/forward.fastqsanger.gz ; [ $? == 0 ] || echo pb $filename/reverse.fastqsanger.gz; done

Vérification du nombre de séquences forward/reverse

Indexation de la référence

Comme nous ne savons pas à quel génome, Salmo salar ? Salmo trutta ? Anguilla rostrata?, ces échantillons sont les plus proches, nous allons générer une heatmap des FASTQ par rapport à ces 3 génomes (et non des échantillons entre eux).

Le génome indexé de Salmo trutta est ici :
/work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/Ref/SalmoTrutta_WholeGenome_GCF_901001165.1_fSalTru1.1_genomic.fna.gz].fasta.gz

Les génomes males et femelles indexés pour Salmo salar sont ici:
/work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/REF/Salmo_salar_Female_AtlanticSalmon_WholeGenome_GCF_000233375.1_ICSASG_v2_genomic.fna.gz].fasta.gz
et
/work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/REF/Salmo_Salar_Male_GCF_905237065.1_Ssal_v3.1_genomic.fna.gz

Le génome indexé pour Anguilla rostrata est ici:
/work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Anguilla_rostrata/Ref/Anguilla_rostrata_2011_GCA_001606085.1_ASM160608v1_genomic.fna].fasta

Aligner sur la référence, filtrer et regrouper

Préparation du fichier 1-mapp.sh pour Salmo sp
Lancement du job:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp $ sarray -c 4 --mem=10G 1-map-SalmoTrutta.sh
Submitted batch job 43989612 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp $ sarray -c 4 --mem=10G 1-map-SalmoSalarMale.sh
Submitted batch job 43989639 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp $ sarray -c 4 --mem=10G  1-map-SalmoSalarFemelle.sh
Submitted batch job 43989666 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp $ sarray -c 4 --mem=10G 1-mapp-AnguillaRostrata.sh
Submitted batch job 43989693 / COMPLETED

Préparation des étapes filter, convert et merge intervals pour Salmo sp

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp $ sbatch 2-filter-convert.sh 
Submitted batch job 43990401 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp $ sbatch -d afterok:43990401 3-merge-intervals.sh 
Submitted batch job 43990402 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp $ sbatch -d afterok:43990402 4-bedtools-coverage.sh 
Submitted batch job 43990403 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp $ sbatch -d afterok:43990403 5-calculate-coverage-by-sample.sh 
Submitted batch job 43990404 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/RESULTATS $ sbatch MATRICE.sh
Submitted batch job 43994757 / COMPLETED
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/RESULTATS $ sed -i 's/; /,/g' MATRICE.txt
(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/RESULTATS $ head MATRICE.txt
Ajouter le header: sample1,sample2,nbnt 
python3 calculate_ratio.py
head MATRICE_ratios.csv avec le header: sample1,sample2,ratio

Code pour générer la heatmap Salmo sp avec les 108 samples

Correspondance entre le nom et le numéro de l'échantillon: Voir le fichier "equivalence_echantillons_Salmo_sp_nb.csv" en bas de cette page wiki.
Salmo trutta = 1 -> 27
Salmo salar male = 28 -> 54
Salmo salar femelle = 55 -> 81
Anguilla rostrata = 82 -> 108

equivalences_Salmo_sp.png

Représentation de la heatmap (aperçu avec seulement une trentaine d'échantillons)

Code pour générer les heatmaps Salmo_sp

Salmo_sp.png

Supprimons les ratios =1 afin de voir si la visualisation est différente pour les échantillons.

$ grep -v ", 1.0" MATRICE_ratios.csv > MATRICE_ratios_without1.csv

Salmo_sp_without1.png

Représentation pour l'ensemble des 108 échantillons

heatmap_all_salmo_sp_without_ratio_1.png

Résultat: Matrice avec l'ensemble des ratios

La matrice est disponible ici: /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/RESULTATS
  • Avec l'ensemble des ratios:MATRICE_ratios.csv
  • Sans les ratios égaux à 1: MATRICE_ratios_without1.csv

6. Compilation des résultats de PALEOMIX

Contact : Stefanie Wagner <>

Voici un lien avec les paths vers la documentation et les résultats (PathToPrefixes: génome de réfrérence utilisés, PathToParameters: paramètre de l'analyses, PathToResults: résultats du mapping + analyse de dégradation):
https://docs.google.com/spreadsheets/d/1jiJQiJA8Ov1tEEs13kdz5g6rINGiTr0b8_cI8x45oic/edit?usp=sharing
Pour la prochaine réunion de travail je propose à l'équipe qui travaille sur la partie criblage en faible profondeur
1/ de compiler les summary files (dans les dossiers de résultats) pour tous les échantillons dans un fichier excel, une ligne par échantillon suivi par les rubriques des summary files :
seq_reads_pairs | seq_trash_pe_1 | seq_trash_pe_1_frac | .... |hits_raw(Anguilla.anguilla-nuc.mt) | ... | hits_length(Salmo.trutta-nuc.mt)
Sample_AUD10996-1
Sample_AUD20212-2
2/ de compiler les pdf de l'analyse de dégradation (Fragmisincorporation_plot.pdf, Length_plot.pdf) pour les données nuc+mito combinés.

Compilation des summary files

Script python pour compiler les summary files

Lancement du script python qui comprends les étapes suivantes: récupération des colonnes Target, Measure, Value de l'ensemble des fichiers .summary des répertoires results_* dans un seul fichier convertit en dataframe. Transposition de la dataframe et suppression des duplicats puis export en fichier csv, avec une ligne par échantillon suivi par les rubriques des summary files.

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX $ python3 compile_results.py

En sortie, le fichier csv est une compilation des résultats des mesures par échantillon: compilation_Resultats.csv Ce fichier est disponible en bas de cette page wiki.
A vérifier: pour faire pivoter la dataframe,j'ai dû enlever les duplicats. Donc vérifier qu'il ne manque pas de données.

compilation_Resultats.png

Les fichiers sont disponibles en bas de page:
compilation_Resultats_2021_mt.csv
compilation_Resultats_2021_comb.csv
compilation_Resultats_2019-2020_mt.csv
compilation_Resultats_2019-2020_comb.csv

Compilation des pdf de l'analyse de dégradation

Question: Compiler les pdf de l'analyse de dégradation (Fragmisincorporation_plot.pdf, Length_plot.pdf) pour les données nuc+mito combinés.

Localisation des fichiers, uniquement dans les répertoires results_20*comb*:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish $ find . -name "Fragmisincorporation_plot.pdf" | head
./results_2021_mt.part2/L5_4_MD4_S36.Homo.sapiens-mt.mapDamage/L5_4_MD4_S36/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L5_12_MD9_S44.Anguilla.rostrata-mt.mapDamage/L5_12_MD9_S44/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L4_1_MB14_S17.Homo.sapiens-mt.mapDamage/L4_1_MB14_S17/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L4_2_MD13_S18.Hucho.hucho-mt.mapDamage/L4_2_MD13_S18/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L5_16_blanc_S48.Homo.sapiens-mt.mapDamage/L5_16_blanc_S48/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L3_9_MB17_S9.Salmo.salar-mt.mapDamage/L3_9_MB17_S9/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L3_2_BC17_S2.Hucho.hucho-mt.mapDamage/L3_2_BC17_S2/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L5_14_BC12_S46.Hucho.hucho-mt.mapDamage/L5_14_BC12_S46/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L4_8_MB4_S24.Homo.sapiens-mt.mapDamage/L4_8_MB4_S24/Fragmisincorporation_plot.pdf
./results_2021_mt.part2/L4_3_MB16_S19.Salmo.trutta-mt.mapDamage/L4_3_MB16_S19/Fragmisincorporation_plot.pdf

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish $ find . -name "Length_plot.pdf" | head
./results_2021_mt.part2/L5_4_MD4_S36.Homo.sapiens-mt.mapDamage/L5_4_MD4_S36/Length_plot.pdf
./results_2021_mt.part2/L5_12_MD9_S44.Anguilla.rostrata-mt.mapDamage/L5_12_MD9_S44/Length_plot.pdf
./results_2021_mt.part2/L4_1_MB14_S17.Homo.sapiens-mt.mapDamage/L4_1_MB14_S17/Length_plot.pdf
./results_2021_mt.part2/L4_2_MD13_S18.Hucho.hucho-mt.mapDamage/L4_2_MD13_S18/Length_plot.pdf
./results_2021_mt.part2/L5_16_blanc_S48.Homo.sapiens-mt.mapDamage/L5_16_blanc_S48/Length_plot.pdf
./results_2021_mt.part2/L3_9_MB17_S9.Salmo.salar-mt.mapDamage/L3_9_MB17_S9/Length_plot.pdf
./results_2021_mt.part2/L3_2_BC17_S2.Hucho.hucho-mt.mapDamage/L3_2_BC17_S2/Length_plot.pdf
./results_2021_mt.part2/L5_14_BC12_S46.Hucho.hucho-mt.mapDamage/L5_14_BC12_S46/Length_plot.pdf
./results_2021_mt.part2/L4_8_MB4_S24.Homo.sapiens-mt.mapDamage/L4_8_MB4_S24/Length_plot.pdf
./results_2021_mt.part2/L4_3_MB16_S19.Salmo.trutta-mt.mapDamage/L4_3_MB16_S19/Length_plot.pdf

Exemples disponibles en bas de page wiki: /work/project/ECOBIOP_ADNa/Data_PaleoFish/results_2021_comb.part1/L10_10_R9_54_e_S122.Anguilla.anguilla-nuc.mt.mapDamage/L10_10_R9_54_e_S122/Fragmisincorporation_plot.pdf Fragmisincorporation_plot.png
et
/work/project/ECOBIOP_ADNa/Data_PaleoFish/results_2021_comb.part1/L10_10_R9_54_e_S122.Anguilla.anguilla-nuc.mt.mapDamage/L10_10_R9_54_e_S122/Length_plot.pdf Length_plot.png

L'information d'intérêt est la substitution du C en T (désamination du C en U puis du U en T), dans le fichier dnacomp.txt:

(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/results_2021_mt.part2/L5_4_MD4_S36.Homo.sapiens-mt.mapDamage/L5_4_MD4_S36 $ more misincorporation.txt
# table produced by mapDamage version 2.0.8
# using mapped file /work/project/ECOBIOP_ADNa/Data_PaleoFish/tmp/70da43fc-33ae-48ff-a80c-821a422e0486/input.
bam and /work/project/ECOBIOP_ADNa/Data_PaleoFish/prefixes.2/Homo.sapiens.mito.fasta as reference file
# Chr: reference from sam/bam header, End: from which termini of DNA sequences, Std: strand of reads
Chr    End    Std    Pos    A    C    G    T    Total    G>A    C>T    A>G    T>C    A>C    A>T    C>G    C>A    T>G    T>A    G>C    G>T    A>-    T>-    C>-    G>-    ->A    ->T    ->C    ->G    S
*    3p    +    1    1    0    0    0    1    0    0    0    0    0    00    0    0    0    0    0    0    0    0    0    0    0    0    0    0
*    3p    +    2    1    0    0    0    1    0    0    0    0    0    00    0    0    0    0    0    0    0    0    0    0    0    0    0    0
*    3p    +    3    0    0    0    1    1    0    0    0    0    0    00    0    0    0    0    0    0    0    0    0    0    0    0    0    0
*    3p    +    4    0    0    1    0    1    0    0    0    0    0    00    0    0    0    0    0    0    0    0    0    0    0    0    0    0
*    3p    +    5    1    0    0    0    1    0    0    0    0    0    00    0    0    0    0    0    0    0    0    0    0    0    0    0    0
*    3p    +    6    0    0    1    0    1    0    0    0    0    0    00    0    0    0    0    0    0    0    0    0    0    0    0    0    0

Script de compilation des résultats: /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX $ more compile_NbSub_RatioSub_by_genome.sh

Lancement du script:

(base) smaman@node104 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX $ ./compile_NbSub_RatioSub_by_genome.sh

Fichier résultat disponible: /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX/compile_nbSub_ratioSub.csv

(base) smaman@node104 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX $ more compile_nbSub_ratioSub.csv
Genome    Sample    NbSub    tot    RatioSub
Anguilla anguilla nuc    Sample_AUD10199-1    76051    11497899    .0066143388
Anguilla rostrata nuc    Sample_AUD10199-1    77357    11425644    .0067704717
Homo sapiens nuc    Sample_AUD10199-1    3022    934759    .0032329188
Hucho hucho nuc    Sample_AUD10199-1    4592    878790    .0052253666
Salmo salar nuc    Sample_AUD10199-1    5125    932003    .0054989093
Salmo trutta nuc    Sample_AUD10199-1    5133    933389    .0054993148
Anguilla anguilla nuc    Sample_AUD10310-2    82383    11587256    .0071097937
Anguilla rostrata nuc    Sample_AUD10310-2    83667    11521411    .0072618709
Homo sapiens nuc    Sample_AUD10310-2    2039    673435    .0030277606
Hucho hucho nuc    Sample_AUD10310-2    3166    605749    .0052265872
...

Vérification :

(base) smaman@node105 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX $ ls RES_compilation_NbSub_RatioSub_by_genome/ | wc -l
2052

(base) smaman@node105 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX $ wc -l compile_nbSub_ratioSub.csv
2053 compile_nbSub_ratioSub.csv

(base) smaman@node105 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX $ head compile_nbSub_ratioSub.csv
Genome    Sample    NbSub    tot    RatioSub
Anguilla anguilla 2019_2020_comb    Sample_AUD10199-1    76051    11497899    .0066143388
Anguilla rostrata 2019_2020_comb    Sample_AUD10199-1    77357    11425644    .0067704717
Homo sapiens 2019_2020_comb    Sample_AUD10199-1    3022    934759    .0032329188
Hucho hucho 2019_2020_comb    Sample_AUD10199-1    4592    878790    .0052253666
Salmo salar 2019_2020_comb    Sample_AUD10199-1    5125    932003    .0054989093
Salmo trutta 2019_2020_comb    Sample_AUD10199-1    5133    933389    .0054993148
Anguilla anguilla 2019_2020_comb    Sample_AUD10310-2    82383    11587256    .0071097937
Anguilla rostrata 2019_2020_comb    Sample_AUD10310-2    83667    11521411    .0072618709
Homo sapiens 2019_2020_comb    Sample_AUD10310-2    2039    673435    .0030277606
...

L’information sur le nombre de substitutions (NbSub) et ce nombre/total substitution (RatioSub) est à chaque fois en fonction d’un génome de référence.

Ce qui signifie qu’il faudrait en sortie :


                            Anguilla anguilla             Anguilla rostrata            Salmo salar                      Salmo trutta                    Homo sapien

                            NbSub  RatioSub            NbSub  RatioSub            NbSub              RatioSub            NbSub  RatioSub            NbSub  RatioSub

N° Echantillon

7. Accès aux résultats

Accès aux matrices contenant les ratios

  • Salmo trutta:
    /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_trutta/RESULTATS/MATRICE_ratios.csv
    
  • Anguilla anguilla:
    /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Anguilla_anguilla/RESULTATS/MATRICE_ratios.csv
    
  • Salmo salar:
    /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_salar/RESULTATS/MATRICE_ratios.csv
    
  • Salmo sp:
    /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/recouvrement_entre_samples_Salmonidae/Salmo_sp/RESULTATS/MATRICE_ratios.csv (limité à 24 échantillons)
    Matrice complète en cours de réalisation.
    

Représentation des heatmaps par période

1/ Pour chaque matrice, remplacer les numéro par le nom de la période + échantillon (en abrégé).
2/ Regrouper toutes la matrices Salmo trutta, Anguilla anguilla, Salmo salar et Salmo sp.
3/ Générer une heatmap pour l'ensemble des données ou par période en fonction de la lisibilité.

Résultats des compilations

Compilation des summary files

Les fichiers sont disponibles en bas de page:
compilation_Resultats_2021_mt.csv
compilation_Resultats_2021_comb.csv
compilation_Resultats_2019-2020_mt.csv
compilation_Resultats_2019-2020_comb.csv

Compilation des PDF


(base) smaman@genologin1 /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX $ head  /work/project/ECOBIOP_ADNa/Data_PaleoFish/SARAH/COMPILATION_PALEOMIX/compile_nbSub_ratioSub.csv
Genome    Sample    NbSub    tot    RatioSub
Anguilla anguilla 2019_2020_comb    Sample_AUD10199-1    76051    11497899    .0066143388
Anguilla rostrata 2019_2020_comb    Sample_AUD10199-1    77357    11425644    .0067704717
Homo sapiens 2019_2020_comb    Sample_AUD10199-1    3022    934759    .0032329188

Est-ce que ce format convient ?

Heatmap par périodes historiques

Message de Natacha:

Voici la table avec les infos sur les ind. dont les périodes historiques : fichier en bas de page Data1-paleofish pour NN-December2022_NameCorrected_PeriodInfoCorrected_WithLizQuinlandData.xlsx
Je me suis trompée avec nos échantillons récents, sur les échantillons anciens nous n’avons pas le sexe puisqu’il s’agit de vertèbres archéologiques

8. Utilisation des primers pour une première catégorisation des échantillons par espèce

En attaché un récapitulatif des primers à utiliser: Amorces_MigrADNe_MV.xlsx, ATL_eels_11_SNP_diagno.xls et Primers_a_utiliser_info.docx

On a les primers universels (MiFish, Tele01, Tele02) mais comme ils sont pas forcément discrimant sur toutes les espèces, il existe les primers Cyp01 (pour les cyprinidés), Gai01 (Gadidés), les primers spécialement développés récemment pour distinguer Anguille, saumon, truite, alose etc. de tous nos poissons d’eau douce en France (c’est la dernière table excel en attaché). Puis les 11 séquences diagnostique pour différencier les 2 espèces anguilles ainsi que les hybrides.

Proposition de pipeline pour rechercher les primers dans les échantillons avec BLAT

1- Reconstruire les séquences avec FLASH en amont de BLAT

FLASH permet de construire des séquences plus grandes en overlappant les paires de bases.

flash -t 8  sample.R1.fastq.gz     sample.R2.fastq.gz     -o sample

Tu obtiens un fastq extended: sample.extendedFrags.fastq sur lequel tu lances ta conversion fastq to fasta puis BLAT de nouveau.

Dans un premier temps, ne lance BLAT que sur un sample (qui est censé contenir un primer) avec un primer court puis avec un primer long.

2- Convertir FASTQ en FASTA avec la commande fastq_to_fasta

#!/bin/bash
#SBATCH -p workq
#SBATCH -t 0-10:60 #Acceptable time formats include "minutes", "minutes:seconds", "hours:minutes:seconds", "days-hours", "days-hours:minutes" and "days-hours:minutes:seconds". 

#Load modules
module load bioinfo/fastx_toolkit-0.0.14

#fastq_to_fasta -i example.fastq -o example.fasta

3- BLAT fasta sample sur primer.fasta

#!/bin/bash
#SBATCH -p unlimitq

#Load modules
module load bioinfo/blatSuite.36

cd /work/project/sigenae/sarah/lncTrout/BLAT/;

#returns no hit
blat -minScore=0 -stepSize=2 database.fa test.fa output_no_hit.psl

#returns 2 hits
blat -minScore=0 -stepSize=1 -t=dna  -q=rna  primer1.fa test.fasta  output_2_hit_primer1.psl

Résultats:
The following information is returned: the score of the alignment, the region of query sequence that matches to the database sequence, the size of the query sequence, the level of identity as a percentage of the alignment and the chromosome and position that the query sequence maps to.

exemple de fichier .psl
match    mis-     rep.     N's    Q gap    Q gap    T gap    T gap    strand    Q            Q       Q        Q      T            T       T        T      block    blockSizes     qStarts     tStarts
         match    match           count    bases    count    bases              name         size    start    end    name         size    start    end    count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
195    5    0    0    0    0    0    0    +    A00924:200:HG5YCDRXY:1:2116:24957:27743    202    1    201    gnl|lncRNA|t178l|29566    277    60    260    1    200,    1,    60,
201    6    0    0    0    0    0    0    +    A00924:200:HG5YCDRXY:2:2242:8368:34147    271    1    208    gnl|lncRNA|t178l|29566    277    61    268    1    207,    1,    61,
202    4    0    0    0    0    0    0    +    A00924:200:HG5YCDRXY:2:2164:11324:19664    256    50    256    gnl|lncRNA|t178l|29566    277    4    210    1    206,    50,    4,
202    4    0    0    0    0    0    0    +    A00924:200:HG5YCDRXY:2:2164:11379:19852    256    50    256    gnl|lncRNA|t178l|29566    277    4    210    1    206,    50,    4,
matches - Number of matching bases that aren't repeats.
misMatches - Number of bases that don't match.
repMatches - Number of matching bases that are part of repeats.
nCount - Number of 'N' bases.
qNumInsert - Number of inserts in query.
qBaseInsert - Number of bases inserted into query.
tNumInsert - Number of inserts in target.
tBaseInsert - Number of bases inserted into target.
strand - defined as + (forward) or - (reverse) for query strand. In mouse, a second '+' or '-' indecates genomic strand.
qName - Query sequence name.
qSize - Query sequence size.
qStart - Alignment start position in query.
qEnd - Alignment end position in query.
tName - Target sequence name.
tSize - Target sequence size.
tStart - Alignment start position in target.
tEnd - Alignment end position in target.
blockCount - Number of blocks in the alignment.
blockSizes - Comma-separated list of sizes of each block.
qStarts - Comma-separated list of start position of each block in query.
tStarts - Comma-separated list of start position of each block in target.

Filtrer avec au moins 99% de couverture, le nombre de fragments retrouvés dans l'échantillon et qui correspondent au primer recherché:

$ cat output_no_hit_SAMPLE_primer1.psl    | awk '{if ((/^[0-9]+/)&&($10!=$14)&&($1/$11>0.99)&&($18==1)) print $10,$14;}'  | wc -l
40

Après, il nous faudra travailer sur ce nombre de fragments pour catégoriser les échantillons par espèce.

Problèmes d'output avec BLAT

Petite description du fichier primers_blat.sh (Disponible parmi les fichiers de ce wiki) pour expliquer ce que j'ai voulu faire et ce qui ne marche pas.

Pour chaque fichier fastq obtenu avec FLASH, stocké dans $qdir, on créé un fichier .fasta.

if file $fastq | grep -q compressed
        then
                gzip -dc $fastq | sed -n '1~4s/^@/>/p;2~p' > $nom.fa
        else
                sed -n '1~4s/^@/>/p;2~4p' $fastq > $nom.fa
        fi

Puis, on BLAT chaque pair forward/reverse de primers (dans un seul fichier $primers) sur ce fichier .fasta et on stock le résultat dans le répertoire $sdir avant de filtrer le nombre de fragments retrouvés (avec ici le seuil fixé à un taux de couverture > 95%).

pri=$(basename -- $primers)
out=$sdir/$nom_$pri_2_hits.psl
blat -minScore=0 -stepSize=1 $primers $nom.fa $out
n=$(cat $out | awk '{if ((/^[0-9]+/)&&($10!=$14)&&($1/$11>0.95)&&($18==1)) print $10,$14;}' | wc -l)
echo "$nom, $pri, $n" 

Le problème, c'est que même si la filtration fonctionne, je ne retrouve pas dans $sdir les fichiers de sortie attendues. Même en essayant de forcer la main au script en écrivant;

blat -minScore=0 -stepSize=1 $primers $nom.fa > $out

Je ne retrouve toujours pas de fichiers .psl dans le répertoire de sauvegarde.

9. Blast+ de nos fichiers nettoyés sur tous les téléostéens de NCBI

Demande: Comme nos échantillons sont des vertèbres archéologiques dont nous ne sommes pas sûr de l’affiliation à l’espèce, on aimerait faire des blasts de nos fichiers nettoyés sur tous les téléostéens de NCBI.

Blast+ est disponible sur genotoul: http://bioinfo.genotoul.fr/index.php/how-to-use/?software=How_to_use_SLURM_NCBI_Blast%2B
TODO: Liste complète des génomes des poissons téléostéens (lien htpp vers la bonne version) pour demander leur mise à disposition et leur indexation sur genologin.

10. Matrice avec les % d’alignement des échantillons sur un génome de référence ?

TODO:
Avec Patrick on se posait la question comment on pouvait récupérer l’information du % d’alignement de notre échantillon sur un génome de référence et la compiler dans un tableau.
Ceci pour les données que tu as fait tourné mais aussi celle que Stéphanie a fait tourné afin de voir les différences.

--> A préciser en réunion ?

paleofish_ratios.png (17.8 KB) Sarah Maman, 10/14/2022 04:40 PM

tableau_correspondance_paleofish_anguilla_anguilla.ods (35.6 KB) Sarah Maman, 10/14/2022 04:47 PM

Salmo_trutta.png (15.3 KB) Sarah Maman, 11/16/2022 02:05 PM

Salmo_salar.png (15.5 KB) Sarah Maman, 11/16/2022 02:10 PM

Salmo_salar_damier.png (15.9 KB) Sarah Maman, 11/16/2022 02:59 PM

Anguilla_anguilla_gauche.png (16.7 KB) Sarah Maman, 11/16/2022 03:09 PM

Salmo_sp.png (14.7 KB) Sarah Maman, 11/21/2022 02:19 PM

Salmo_sp_without1.png (17.2 KB) Sarah Maman, 11/21/2022 02:25 PM

compilation_Resultats.png (160 KB) Sarah Maman, 11/23/2022 11:00 AM

Length_plot.pdf (11 KB) Sarah Maman, 11/23/2022 11:22 AM

Fragmisincorporation_plot.pdf (19.3 KB) Sarah Maman, 11/23/2022 11:22 AM

Fragmisincorporation_plot.png (85.2 KB) Sarah Maman, 11/23/2022 11:24 AM

Length_plot.png (97.3 KB) Sarah Maman, 11/23/2022 11:24 AM

avancement_novembre2022.odp (821 KB) Sarah Maman, 11/23/2022 11:46 AM

compilation_Resultats_2021_mt.csv Magnifier (74.7 KB) Sarah Maman, 11/23/2022 04:16 PM

compilation_Resultats_2021_comb.csv Magnifier (90.5 KB) Sarah Maman, 11/23/2022 04:16 PM

compilation_Resultats_2019-2020_mt.csv Magnifier (20.4 KB) Sarah Maman, 11/23/2022 04:16 PM

compilation_Resultats_2019-2020_comb.csv Magnifier (24.7 KB) Sarah Maman, 11/23/2022 04:16 PM

Amorces_MigrADNe_MV.xlsx (35.8 KB) Sarah Maman, 01/17/2023 04:06 PM

ATL_eels_11_SNP_diagno.xls (51 KB) Sarah Maman, 01/17/2023 04:06 PM

Primers_a_utiliser_info.docx (19 KB) Sarah Maman, 01/17/2023 04:06 PM

Data1-paleofish pour NN-December2022_NameCorrected_PeriodInfoCorrected_WithLizQuinlandData.xlsx (597 KB) Sarah Maman, 01/18/2023 04:06 PM

equivalence_echantillons_Salmo_sp_nb.csv Magnifier (4.03 KB) Sarah Maman, 01/20/2023 09:53 AM

equivalences_Salmo_sp.png (19.8 KB) Sarah Maman, 01/20/2023 09:53 AM

heatmap_all_salmo_sp_without_ratio_1.png (28.5 KB) Sarah Maman, 01/20/2023 04:21 PM

primers_blat.sh Magnifier - script blat primers teleostes janvier (1.73 KB) Patrick Jacques, 01/31/2023 02:17 PM