Wiki¶
- Wiki
- CR réunions
- 1. Description du projet
- 2. Pipeline pour Anguilla anguilla
- 3. Pipeline pour le saumon Salmonidae Salmo salar
- 4. Pipeline pour le saumon Salmonidae Salmo trutta
- 5. Pipeline pour le saumon Salmonidae Salmo sp
- 6. Compilation des résultats de PALEOMIX
- 7. Accès aux résultats
- 8. Utilisation des primers pour une première catégorisation des échantillons par espèce
- 9. Blast+ de nos fichiers nettoyés sur tous les téléostéens de NCBI
- 10. Matrice avec les % d’alignement des échantillons sur un génome de référence ?
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 <natacha.nikolic@inrae.fr>
Stefanie Wagner <stefwag2002@gmail.com>
Aurelie Manicki <aurelie.manicki@inrae.fr>
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¶
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
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
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 |
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.
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.0Correspondance 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
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
Représentation de la heatmap (aperçu avec seulement une trentaine d'échantillons)¶
Code pour générer les heatmaps Salmo_sp
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
Représentation pour l'ensemble des 108 échantillons¶
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 <stefanie.wagner@inrae.fr>
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.
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
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
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)
tableau_correspondance_paleofish_anguilla_anguilla.ods (35.6 KB)
Salmo_trutta.png (15.3 KB)
Salmo_salar.png (15.5 KB)
Salmo_salar_damier.png (15.9 KB)
Anguilla_anguilla_gauche.png (16.7 KB)
Salmo_sp.png (14.7 KB)
Salmo_sp_without1.png (17.2 KB)
compilation_Resultats.png (160 KB)
Length_plot.pdf (11 KB)
Fragmisincorporation_plot.pdf (19.3 KB)
Fragmisincorporation_plot.png (85.2 KB)
Length_plot.png (97.3 KB)
avancement_novembre2022.odp (821 KB)
compilation_Resultats_2021_mt.csv 
(74.7 KB)
compilation_Resultats_2021_comb.csv
(90.5 KB)
compilation_Resultats_2019-2020_mt.csv
(20.4 KB)
compilation_Resultats_2019-2020_comb.csv
(24.7 KB)
Amorces_MigrADNe_MV.xlsx (35.8 KB)
ATL_eels_11_SNP_diagno.xls (51 KB)
Primers_a_utiliser_info.docx (19 KB)
Data1-paleofish pour NN-December2022_NameCorrected_PeriodInfoCorrected_WithLizQuinlandData.xlsx (597 KB)
equivalence_echantillons_Salmo_sp_nb.csv
(4.03 KB)
equivalences_Salmo_sp.png (19.8 KB)
heatmap_all_salmo_sp_without_ratio_1.png (28.5 KB)
primers_blat.sh
- script blat primers teleostes janvier
(1.73 KB)