|
|
**Stratégie** : cutadapt, sickle (optionnel), et bwa sur le génome humain, on ne garde que les paires de reads ne mappant pas.
|
|
|
|
|
|
# Nettoyage
|
|
|
|
|
|
## cutadapt
|
|
|
|
|
|
cutadapt version ???
|
|
|
```
|
|
|
PATH/cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -o Sample_R1.trimmed.fastq -p Sample_R2.trimmed.fastq PATH/RawData/Sample_R1.fastq.gz PATH/RawData/Sample_R2.fastq.gz
|
|
|
```
|
|
|
|
|
|
# sickle
|
|
|
|
|
|
sickle --version
|
|
|
sickle version 1.210
|
|
|
|
|
|
```
|
|
|
sickle pe -f Sample_R1.trimmed.fastq -r Sample_R2.trimmed.fastq -t sanger -o Sample_R1.trimmedQual.fastq -p Sample_R2.trimmedQual.fastq -s Sample_Single.trimmedQual.fastq
|
|
|
```
|
|
|
|
|
|
|
|
|
# alignement
|
|
|
|
|
|
bwa humain (hg19) des reads trimmées pour la qualité : à mettre en option pas toujours nécessaire**
|
|
|
|
|
|
```
|
|
|
PATH/bwa mem PATH_BANK/ensembl_homo_sapiens_genome Sample_R1.trimmedQual.fastq Sample_R2.trimmedQual.fastq -t 4 | samtools view -bhS - > Sample.bam
|
|
|
```
|
|
|
|
|
|
Lancement cluster
|
|
|
```
|
|
|
qarray -pe parallel_smp 4 -l mem=30G -l h_vmem=33G bwa.sh
|
|
|
samtools flagstat Sample.bam > Sample.stat
|
|
|
```
|
|
|
|
|
|
## Filtre
|
|
|
|
|
|
**je ne garde que les paires non mappées**
|
|
|
|
|
|
```
|
|
|
PATH/samtools view -b -f 12 Sample.bam > SampleWithoutPairedUnmapped.bam
|
|
|
```
|
|
|
|
|
|
```
|
|
|
PATH/samtools flagstat SampleWithoutPairedUnmapped.bam > SampleWithoutPairedUnmapped.stat
|
|
|
```
|
|
|
|
|
|
**je reconstruis les fastq**
|
|
|
|
|
|
```
|
|
|
PATH/bamToFastq -i SampleWithoutPairedUnmapped.bam -fq SampleWithoutHumanR1.fastq -fq2 SampleWithoutHumanR2.fastq
|
|
|
```
|
|
|
|
|
|
**compresser**
|
|
|
```
|
|
|
gzip SampleWithoutHumanR1.fastq
|
|
|
gzip SampleWithoutHumanR2.fastq
|
|
|
```
|
|
|
|
|
|
**fastQC**
|
|
|
```
|
|
|
PATH/fastqc -o FastqcDir --casava --nogroup Sample1WithoutHumanR1.fastq.gz Sample2WithoutHumanR1.fastq.gz Sample1WithoutHumanR2.fastq.gz Sample2WithoutHumanR1.fastq.gz Sample3WithoutHumanR2.fastq.gz Sample3WithoutHumanR1.fastq.gz Sample4WithoutHumanR2.fastq.gz Sample4WithoutHumanR2.fastq.gz
|
|
|
```
|
|
|
|
|
|
########################################################################################################
|
|
|
|
|
|
**En ce qui concerne les options ou informations à sortir du composant**
|
|
|
|
|
|
Je pense au nom des échantillons dans notre cas : Sample par exemple.
|
|
|
|
|
|
Il faudrait permettre de choisir s’il y a nécessité d’enlever les contaminants (ici les reads humains) et de spécifier le génome du contaminant (ici hg19)
|
|
|
|
|
|
########################################################################################################
|
|
|
|
|
|
# Assemblage et annotation PROKKA
|
|
|
|
|
|
## 1/ Metaspades
|
|
|
|
|
|
### a/ Lancements individuels des lots
|
|
|
```
|
|
|
for i in samples
|
|
|
do
|
|
|
echo "PATH/metaspades.py -1 "$i"WithoutHumanR1.fastq.gz -2 "$i"WithoutHumanR2.fastq.gz -o "$i"_metaspades" >> metaspades.sh
|
|
|
done
|
|
|
```
|
|
|
|
|
|
**Traitement en parallèle**
|
|
|
```
|
|
|
qarray -l mem=4G,h_vmem=256G -q smpq metaspades.sh
|
|
|
```
|
|
|
|
|
|
## 2/ Megahit
|
|
|
|
|
|
### a/ Lancements individuels des lots
|
|
|
|
|
|
**Contenu du fichier megahit.sh**
|
|
|
```
|
|
|
PATH/megahit -t 48 -1 SampleWithoutHumanR1.fastq.gz -2 SampleWithoutHumanR2.fastq.gz -o Sample
|
|
|
```
|
|
|
|
|
|
**Traitement en série**
|
|
|
```
|
|
|
qsub -q smpq -l mem=24G,h_vmem=512G megahit.sh
|
|
|
```
|
|
|
## 3/ Annotation prokka
|
|
|
|
|
|
**Exemple pour un lot**
|
|
|
```
|
|
|
PATH/prokka --outdir PROKKA_Sammple --prefix Sample Sample_metaspades_scaffolds.fasta --centre X –compliant
|
|
|
```
|
|
|
|
|
|
## 4/ CD-HIT individuels
|
|
|
|
|
|
### Clustering
|
|
|
|
|
|
Le clustering est fait sur les CDS en sortie de prokka pour éliminer les CDS redondantes (à 5% = même espèce)
|
|
|
```
|
|
|
cd-hit-est -c 0.95 -i Sample.ffn -o Sample.faa.cd-hit-est.095
|
|
|
```
|
|
|
|
|
|
**NB.** Clustering chois à 95% sur les CDS en nucléotides.
|
|
|
|
|
|
**NB2.** le fichier clstr (cluster) donne le lien entre le représentant du cluster et ses membres. Le script python ci-dessous permet de récuperer un fichier tabulé contenant cette information de lien.
|
|
|
|
|
|
|
|
|
### Script python (cd-hit_produce_table_clstr.py) pour récupérer les membres de chaque cluster
|
|
|
```
|
|
|
cat fic.clstr | python cd-hit_produce_table_clstr.py | more
|
|
|
|
|
|
#!/bin/env python
|
|
|
import sys, re
|
|
|
|
|
|
#init dictionaries:
|
|
|
ref = ""
|
|
|
seqs = []
|
|
|
|
|
|
while 1 :
|
|
|
line = sys.stdin.readline()
|
|
|
#print line
|
|
|
if line == '' :
|
|
|
break
|
|
|
else :
|
|
|
if line[0] == ">":
|
|
|
for seq in seqs :
|
|
|
print ref+"\t"+seq
|
|
|
ref = ""
|
|
|
seqs = []
|
|
|
else:
|
|
|
a, b = line.split('>', 1)
|
|
|
name = b.split("...")[0]
|
|
|
rep = (line.rstrip()[-1] == '*')
|
|
|
if rep :
|
|
|
ref = name
|
|
|
seqs.append(name)
|
|
|
else :
|
|
|
seqs.append(name)
|
|
|
|
|
|
for seq in seqs :
|
|
|
print ref+"\t"+seq
|
|
|
```
|
|
|
|
|
|
## 5/ CD-HIT global (synthèse de tous les lots)
|
|
|
|
|
|
**Concatenation de tous les fasta des cd-hit individuels avec le même paramètre**
|
|
|
```
|
|
|
cat *cd-hit-est.095 > All-cd-hit-est.095
|
|
|
cd-hit-est -c 0.95 -i All-cd-hit-est.095 -o All-cd-hit_cd-hit-est.095
|
|
|
```
|
|
|
|
|
|
## 6/ Récupération pour chaque protéine initiale du cluster dans lequel elle figure
|
|
|
|
|
|
Faire la jointure des tables pour avoir le lien entre chaque protéines initiales et les protéines finales. Pour faire un fichier de référence des protéines pour l'annotation.
|
|
|
|
|
|
|
|
|
## 7/ Quantification des reads alignés de chaque gène pour chaque condition
|
|
|
|
|
|
**2 solutions** : soit aligner les séquences sur les gènes soit les aligner sur les contigs et utiliser le gff (transformé en gtf) pour extraire les quantifications par gène avec featureCounts.
|
|
|
|
|
|
**Pour aligner sur les gènes**
|
|
|
- récupération des gènes correspondant aux protéines du set final (utilisation du script ~sigenae/bin/fasta_extract.pl d'un fichier de tous les transcrits)
|
|
|
- indexation du fichier avec bwa index
|
|
|
- alignement des séquences des différents lots avec bwa mem
|
|
|
- transformer le fichier sam en bam avec samtools view et samtools sort et samtools index
|
|
|
- récupération des quantifications avec samtools idxstats
|
|
|
|
|
|
**Pour aligner sur les contigs**
|
|
|
- transformation du gff en gtf avec gffread
|
|
|
```
|
|
|
git clone https://github.com/gpertea/gclib
|
|
|
git clone https://github.com/gpertea/gffread
|
|
|
cd gffread/
|
|
|
make release
|
|
|
```
|
|
|
- indexation du fichier des contigs avec bwa index
|
|
|
- alignement des séquences des différents lots avec bwa mem
|
|
|
- transformer le fichier sam en bam avec samtools view et samtools sort et samtools index
|
|
|
- quantification des transcrits avec featureCounts
|
|
|
|
|
|
**Exemple pour featureCounts :**
|
|
|
```
|
|
|
featureCounts -s 1 -O -p -t exon -g gene_id -a sample1.gtf -o sample1.featureCount.count sample1.sort.bam
|
|
|
```
|
|
|
## 8/ Synthèse des quantifications
|
|
|
|
|
|
A partir des fichiers de quantification des différents conditions et du fichier de lien entre les transcrits du set final et les transcrits de chaque condition, refaire un tableau dans lequel chaque ligne corresponde à un transcrit du set final, chaque colonne correspond à une condition et à la croisée des lignes et des colonne on trouve la somme des nombre de séquences des transcrits qui correspondent pour cette condition au transcrit du set de référence.
|
|
|
|
|
|
**J'avais fait ça en utilisant des commandes shell, mais il vaudrait mieux faire un programme python qui produise le tableau.**
|
|
|
|
|
|
########################################################################################################
|
|
|
|
|
|
# Classification taxonomique avec Kaiju
|
|
|
|
|
|
- Versions des logiciels utilisés : Kaiju 1.6.3 (makeDB.sh, kaiju, kaiju2krona, kaijuReport, mergeOutputs), KronaTools-2.7 (ktImportText), R-3.5.1 (histogrammes).
|
|
|
|
|
|
- Version de la database utilisée pour kaiju (makeDB.sh) : Complete Reference Genomes from NCBI RefSeq (20 février 2019)
|
|
|
|
|
|
## Downloading and compiling Kaiju 1.6.3
|
|
|
|
|
|
```
|
|
|
git clone https://github.com/bioinformatics-centre/kaiju.git
|
|
|
cd kaiju/src
|
|
|
make
|
|
|
PATH="/work/jfourquet/antiselfish/data_Analyse_Kaiju/kaiju/bin/:$PATH"
|
|
|
mkdir kaijudb
|
|
|
cd kaijudb
|
|
|
```
|
|
|
|
|
|
## Création de la database et de l'index
|
|
|
|
|
|
### Dans `script_index_kaiju.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 24:00:00
|
|
|
#SBATCH --cpus-per-task=5
|
|
|
#SBATCH --mem=40GB
|
|
|
|
|
|
makeDB.sh -r -t 5
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_index_kaiju.sh`
|
|
|
|
|
|
## Kaiju MEM
|
|
|
|
|
|
### Dans `script_kaiju_MEM.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 24:00:00
|
|
|
#SBATCH --cpus-per-task=25
|
|
|
|
|
|
kaiju -z 25 -t nodes.dmp -f kaiju_db.fmi -i ../SampleWithoutHumanR1.fastq.gz -j ../SampleWithoutHumanR2.fastq.gz -o SampleWithoutHuman_kaiju_MEM.out -a mem
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_kaiju_MEM.sh`
|
|
|
|
|
|
## Kaiju MEM : mode verbose
|
|
|
**Pour avoir la longueur maximale de correspondance entre le read et un élément de la database**
|
|
|
|
|
|
### Dans `script_kaiju_MEM_option_v.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 10:00:00
|
|
|
#SBATCH --cpus-per-task=25
|
|
|
|
|
|
kaiju -z 25 -t nodes.dmp -f kaiju_db.fmi -i ../SampleWithoutHumanR1.fastq.gz -j ../SampleWithoutHumanR2.fastq.gz -o SampleWithoutHuman_kaiju_MEM_verbose.out -a mem -v
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_kaiju_MEM_option_v.sh`
|
|
|
|
|
|
## Kaiju MEM : génération kronas
|
|
|
|
|
|
### Formattage des fichiers kaiju MEM pour pouvoir générer des kronas
|
|
|
|
|
|
#### Dans `script_krona_kaiju_MEM_option_u.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 1:00:00
|
|
|
|
|
|
kaiju2krona -t nodes.dmp -n names.dmp -i SampleWithoutHuman_kaiju_MEM.out -o SampleWithoutHuman_kaiju_MEM_without_unassigned.out.krona -u
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
sbatch `script_krona_kaiju_MEM_option_u.sh`
|
|
|
|
|
|
### Génération des kronas
|
|
|
|
|
|
#### Dans `script_krona_without_unassigned_MEM.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 1:00:00
|
|
|
|
|
|
module purge
|
|
|
|
|
|
module load bioinfo/KronaTools-2.7
|
|
|
|
|
|
ktImportText -o SampleWithoutHuman_kaiju_MEM_without_unassigned.out.krona.html SampleWithoutHuman_kaiju_MEM_without_unassigned.out.krona
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_krona_without_unassigned_MEM.sh`
|
|
|
|
|
|
|
|
|
## Kaiju MEM : lancement de kaijuReport
|
|
|
**Pour avoir un summary, ici exemple pour phylum mais peut être généré aussi pour class, order, family, genus et species.**
|
|
|
|
|
|
### Dans `script_kaiju_summary_phylum_MEM.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 1:00:00
|
|
|
|
|
|
kaijuReport -t nodes.dmp -n names.dmp -i SampleWithoutHuman_kaiju_MEM.out -r phylum -o SampleWithoutHuman_kaiju_MEM.out.phylum_summary
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_kaiju_summary_phylum_MEM.sh`
|
|
|
|
|
|
|
|
|
## Kaiju MEM : merge avec fichiers kraken
|
|
|
|
|
|
### Sort fichiers kaiju
|
|
|
**Obligatoire avant merge avec fichiers kraken**
|
|
|
|
|
|
#### Dans `Sort_kaijuMEM_files.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 24:00:00
|
|
|
|
|
|
sort -k2,2 SampleWithoutHuman_kaiju_MEM.out >SampleWithoutHuman_kaiju_MEM.out.sort
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch Sort_kaijuMEM_files.sh`
|
|
|
|
|
|
### Sort fichiers kraken
|
|
|
**Obligatoire avant merge avec fichiers kaiju**
|
|
|
|
|
|
#### Dans `Sort_kraken_files.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 24:00:00
|
|
|
|
|
|
sort -k2,2 krakenReadsSample_Filtered0.1 >krakenReadsSample_Filtered0.1.sort
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch Sort_kraken_files.sh`
|
|
|
|
|
|
### Merge fichiers kaiju et kraken triés (avec valeurs par défaut des options)
|
|
|
|
|
|
#### Dans `Kaiju_merge_kraken_noOption.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 1:00:00
|
|
|
|
|
|
mergeOutputs -i SampleWithoutHuman_kaiju_MEM.out.sort -j krakenReadsSample_Filtered0.1.sort -o Sample_combined_noOption.out -v
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch Kaiju_merge_kraken_noOption.sh`
|
|
|
|
|
|
### Diagrammes de Venn sous R
|
|
|
```
|
|
|
library(VennDiagram)
|
|
|
|
|
|
# Sample
|
|
|
grid.newpage()
|
|
|
venn.plot <- draw.pairwise.venn(49379485, 9101807, 8229754, c("Kaiju", "Kraken"), fill = c("#DE2916", "#318CE7"), col = c("#DE2916", "#318CE7"), alpha=0.75, lwd=c(1,1), cat.cex=c(2.5,2.5), cat.pos=c(0,90), cex=c(2.5,2.5,2.5), margin=0.05, cat.dist=c(0.02,0.04), ext.pos=25, ext.length=0.75);
|
|
|
```
|
|
|
|
|
|
## Kaiju MEM : obtention des histogrammes de la distribution des reads classifiés selon la longueur maximale du read trouvée par kaiju MEM
|
|
|
|
|
|
#### Dans `Script_histogramme.R` :
|
|
|
|
|
|
```
|
|
|
#!/usr/bin/env Rscript
|
|
|
args = commandArgs(trailingOnly=TRUE)
|
|
|
|
|
|
# test if there is at least one argument: if not, return an error
|
|
|
if (length(args)==0) {
|
|
|
stop("At least one argument must be supplied (input file).n", call.=FALSE)
|
|
|
} else if (length(args)>=1) {
|
|
|
for (i in 1:length(args))
|
|
|
{
|
|
|
tab_init <- read.table(args[i], header=FALSE, sep="\t", fill=TRUE, col.names = paste0("V",seq_len(7)))
|
|
|
print(head(tab_init))
|
|
|
tab <- tab_init[tab_init$V1=="C",]
|
|
|
print(head(tab))
|
|
|
tab_hist <- hist(tab$V4, breaks=seq(11,50,1), plot=FALSE)
|
|
|
print(paste0(args[i], " breaks and counts sum"))
|
|
|
print(tab_hist$breaks)
|
|
|
print(sum(tab_hist$counts))
|
|
|
jpeg(paste0(args[i], '_Kaiju_MEM_counts.jpg'))
|
|
|
plot(tab_hist, ylim=c(0,12000000))
|
|
|
dev.off()
|
|
|
tab_hist$counts <- tab_hist$counts/sum(tab_hist$counts)
|
|
|
jpeg(paste0(args[i], "_Kaiju_MEM_v2.jpg"))
|
|
|
plot(tab_hist, ylim=c(0,0.25))
|
|
|
dev.off()
|
|
|
print(paste0(args[i], " breaks and counts sum"))
|
|
|
print(tab_hist$breaks)
|
|
|
print(sum(tab_hist$counts))
|
|
|
}
|
|
|
}
|
|
|
|
|
|
```
|
|
|
|
|
|
### Lancement script
|
|
|
|
|
|
```
|
|
|
srun --pty bash
|
|
|
|
|
|
module purge
|
|
|
|
|
|
module load system/R-3.5.1
|
|
|
|
|
|
Rscript Script_histogramme.R SampleWithoutHuman_kaiju_MEM_verbose.out
|
|
|
```
|
|
|
|
|
|
|
|
|
## Kaiju Greedy (5 mismatchs)
|
|
|
|
|
|
### Dans `script_kaiju_greedy.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 24:00:00
|
|
|
#SBATCH --cpus-per-task=25
|
|
|
|
|
|
kaiju -z 25 -t nodes.dmp -f kaiju_db.fmi -i SampleWithoutHumanR1.fastq.gz -j SampleWithoutHumanR2.fastq.gz -a greedy -e 5 -o SampleWithoutHuman_kaiju_greedy_e5.out
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_kaiju_greedy.sh`
|
|
|
|
|
|
## Kaiju Greedy (5 mismatchs) : summary
|
|
|
**Ici exemple pour phylum mais peut être généré aussi pour class, order, family, genus et species.**
|
|
|
|
|
|
### Dans `script_kaiju_summary_phylum_greedy_e5.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 1:00:00
|
|
|
|
|
|
kaijuReport -t nodes.dmp -n names.dmp -i SampleWithoutHuman_kaiju_greedy_e5.out -r phylum -o SampleWithoutHuman_kaiju_greedy_e5.out.phylum_summary
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_kaiju_summary_phylum_greedy_e5.sh`
|
|
|
|
|
|
## Kaiju Greedy (5 mismatchs) : génération du krona
|
|
|
|
|
|
### Formattage du fichier pour génération du krona
|
|
|
|
|
|
#### Dans `script_krona_kaiju_greedy_e5_option_u.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 1:00:00
|
|
|
|
|
|
kaiju2krona -t nodes.dmp -n names.dmp -i SampleWithoutHuman_kaiju_greedy_e5.out -o SampleWithoutHuman_kaiju_greedy_e5_without_unassigned.out.krona -u
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_krona_kaiju_greedy_e5_option_u.sh`
|
|
|
|
|
|
### Génération du krona
|
|
|
|
|
|
#### Dans `script_krona_greedy_e5_without_unassigned.sh` :
|
|
|
|
|
|
```
|
|
|
#!/bin/bash
|
|
|
#SBATCH -p workq
|
|
|
#SBATCH -t 1:00:00
|
|
|
|
|
|
module purge
|
|
|
|
|
|
module load bioinfo/KronaTools-2.7
|
|
|
|
|
|
ktImportText -o SampleWithoutHuman_kaiju_greedy_e5_without_unassigned.out.krona.html SampleWithoutHuman_kaiju_greedy_e5_without_unassigned.out.krona
|
|
|
```
|
|
|
Lancement sur slurm :
|
|
|
`sbatch script_krona_greedy_e5_without_unassigned.sh` |
|
|
\ No newline at end of file |