Protocole analyse données NGS » History » Version 6

« Previous - Version 6/8 (diff) - Next » - Current version
Christine Tranchant, 03/02/2015 05:59 PM


Protocole analyse données NGS (Christine)

NGS DNA analysis

1- Quelques stats

Pour avoir la taille des fichiers

du -sh *fastq 
311M    C3KB2ACXX_5_12_11_debar.fastq
525M    C3KB2ACXX_5_12_12_debar.fastq

Pour avoir le nombre de séquences

wc -l *.fastq | awk '{ print $2" \t "$1/4}'
C3KB2ACXX_5_12_11_debar.fastq      1354891
C3KB2ACXX_5_12_12_debar.fastq      2249353
C3KB2ACXX_5_12_1_debar.fastq      1231436
C3KB2ACXX_5_12_3_debar.fastq      953312
C3KB2ACXX_5_12_7_debar.fastq      3144855
C3KB2ACXX_5_13_2_debar.fastq      0

quelques stats avec ea-utils

https://code.google.com/p/ea-utils/

fastq-stats -D /data/projects/coffea-can/GBS/ftp.solgenomics.net/mixed_samples/GBS/fastq/C3KB2ACXX_7_10_11_debar.fastq
reads    3139441
len    93
len mean    93.0000
len stdev    0.0000
len min    93
phred    33
window-size    2000000
cycle-max    35
qual min    2
qual max    41
qual mean    35.8426
qual stdev    9.4216
%A    24.5623
%C    26.8162
%G    23.9754
%T    24.4093
%N    0.2367
total bases    291968013

Plus long avec l'option -D et -x pour avoir des stats par position

[tranchant@node12 ~]$ fastq-stats /data/projects/coffea-can/GBS/ftp.solgenomics.net/mixed_samples/GBS/fastq/C3KB2ACXX_7_10_11_debar.fastq 
reads    3139441
len    93
len mean    93.0000
len stdev    0.0000
len min    93
phred    33
window-size    2000000
cycle-max    35
dups    1839623
%dup    58.5972
unique-dup seq    235189
min dup count    2
dup seq     1    13508    CAGCACGACACTTGTATTGCTCTCCCACAACCCCG
dup seq     2    13106    CAGCCTAAACCGTGAAAACGGGGTTGTGGGAGAGC
dup seq     3    7948    CTGCGTTCGGGAAGGATGAATCGCTCCCGAAAAGG
dup seq     4    7884    CTGCGTCTCTCCCGAGGCCGGGACCTTTTATCATG
dup seq     5    7839    CAGCATTCCGAGTAACTCCTCAACCGGGAGTTCCA
dup seq     6    7618    CAGCTACCGCGGCCCCCGCTTCTTCAGGTGGAACT
dup seq     7    7508    CTGCTGGATTGGCTGTAGGGCTTGCTTCTATTGGA
dup seq     8    7055    CTGCTGGCACAGAGTTAGCCGATGCTTATTCCCCA
dup seq     9    5667    CTGCTGAGGGTAATGAAATTATCCGCGAGGCTAGT
dup seq     10    5577    CAGCGCCCAATACACCGGCAACTCCCATCATATGA
dup mean    8.8219
dup stddev    75.7727
qual min    2
qual max    41
qual mean    35.8426
qual stdev    9.4216
%A    24.5623
%C    26.8162
%G    23.9754
%T    24.4093
%N    0.2367
total bases    291968013

2 - fastqc

http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
  • fastqc *.fastq -o </tt>.

3- cutadapt

https://code.google.com/p/cutadapt/
revseq obligatary
  • cutadapt $(<cutadapt.conf) file.fastq -o file.CUTADAPT.fastq
    
    /usr/local/cutadapt-1.2.1/bin/cutadapt  -b AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT -b ... -q 20 -O 10 -m 35  file.fastq -o file.CUTADAPT.fastq

4- repairing

  • perl pairing.pl file_1.CUTADAPT.fastq file_2.CUTADAPT.fastq file_1.REPAIRING.fastq file_2.REPAIRING.fastq file_single.REPAIRING.fastq
    

5- bwa

http://bio-bwa.sourceforge.net/

bwa index

  • bwa index -a is

Rk : is for short genome and bwtsw is for genome >2Gb

bwa aln

  • bwa aln reference file_forward.fastq > file_forward.sai
    
    bwa aln reference file_reverse.fastq > file_reverse.sai
    

bwa sampe / samse

  • bwa sampe reference file_forward.sai file_reverse.sai  file_forward.fastq file_reverse.fastq -f file.sam  -r '@RG        ID:RG  SM:RG  PL:Illumina'
    
    bwa samse reference file.sai file.fastq file_reverse.fastq -f file.sam  -r '@RG        ID:RG  SM:RG  PL:Illumina'

6- picardtools

http://broadinstitute.github.io/picard/

CreateSequenceDictionary

  • /usr/local/java/latest/bin/java -Xmx8g  -jar /usr/local/picard-tools-1.83/CreateSequenceDictionary.jar REFERENCE=../Reference.fasta OUTPUT=../Reference.dict

SortSam

  • /usr/local/java/latest/bin/java -Xmx12g -jar /home/sabotf/sources/picard-tools/SortSam.jar SORT_ORDER=coordinate VALIDATION_STRINGENCY=SILENT INPUT=file.BWASAMPE.sam OUTPUT=file.PICARDTOOLSSORT.bam
    

Do the createSequenceDictionnary on the reference before

7- samtools

http://samtools.sourceforge.net/

samtools index

  • samtools index file.bam

samtools view

  • samtools view -h -b -f=0*02 -o file.SAMTOOLSVIEW-MAPPED.bam file.PICARDTOOLSSORT.bam
  • samtools view -h -b -F=0*02 -o file.SAMTOOLSVIEW-UNMAPPED.bam file.PICARDTOOLSSORT.bam
  • samtools view -h -b -F=4 -o file.SAMTOOLSVIEW-MAPPED.bam file.PICARDTOOLSSORT.bam

samtools index

samtools faidx

  • samtools faidx reference.fasta

8- picardToold reorder

http://broadinstitute.github.io/picard/command-line-overview.html#ReorderSam

  • /usr/local/java/latest/bin/java -Xmx12g -jar /home/sabotf/sources/picard-tools/ReorderSam.jar INPUT=20135.SAMTOOLSVIEW-REHEADER.bam OUTPUT=20135.PICARDTOOLSREORDER.bam REFERENCE=Coffea_pseudomolecules.fa
    

Rk: indispensable if you have the error message with gatk realignerTargetCreator
input file reads and reference have incompatible contigs

9- picardtools AddOrReplaceReadGroups

  • /usr/local/java/latest/bin/java -Xmx12g -jar /home/sabotf/sources/picard-tools/AddOrReplaceReadGroups.jar INPUT=20135.SAMTOOLSVIEW-MAPPED.bam OUTPUT=20135.PICARDTOOLSADDORREPLACEREADGROUPS.bam RGPL=Illumina RGLB=20135 RGPU=20135 RGSM=20135 CREATE_INDEX=True VALIDATION_STRINGENCY=LENIENT > 20135.PICARDTOOLSADDORREPLACEREADGROUPS.log
    

RGPL= Read Group Platform (i.e. Illumina)
RGLB= Read Group Library (i.e. DNA preparation library identity i.e "HG19")
RGPU= Read Group Platform Unit (i.e Barcode)
RGSM= Read Group Sample Name (i.e. Sample Name)

Rk: ?? picard error: Exception in thread "main" net.sf.samtools.SAMFormatException: SAM validation error: ERROR: Read name CIGAR M operator maps off end of reference ??
http://seqanswers.com/forums/archive/index.php/t-16905.html
Use option => VALIDATION_STRINGENCY=LENIENT

Rk: indispensable for running gatk realignerTargetCreator

10- samtools index

11- gatk

https://www.broadinstitute.org/gatk/

gatk realignerTargetCreator

  • /usr/local/java/latest/bin/java -Xmx12g -jar /home/sabotf/sources/GenomeAnalysisTK/GenomeAnalysisTK.jar -T RealignerTargetCreator --fix_misencoded_quality_scores -fixMisencodedQuals -R ../Reference.fasta -I ../RC3.SAMTOOLSVIEW.bam -o ../RC3.GATKREALIGNERTARGETCREATOR.intervals
    
    /usr/java/jre1.7.0_51/bin/java -Xmx12g -jar /usr/local/GenomeAnalysisTK-3.3/GenomeAnalysisTK.jar -T RealignerTargetCreator -nt 4  -R ../Reference.fasta -I ../RC3.SAMTOOLSVIEW.bam -o ../RC3.GATKREALIGNERTARGETCREATOR.intervals
    

gatk indelRealigner

  • /usr/local/java/latest/bin/java -Xmx12g -jar /home/sabotf/sources/GenomeAnalysisTK/GenomeAnalysisTK.jar -T IndelRealigner --fix_misencoded_quality_scores -fixMisencodedQuals -R ../Reference.fasta -I ../RC3.SAMTOOLSVIEW.bam -targetIntervals ../RC3.GATKREALIGNERTARGETCREATOR.intervals -o ../RC3.GATKINDELREALIGNER.bam
    
    /usr/java/jre1.7.0_51/bin/java -Xmx12g -jar /usr/local/GenomeAnalysisTK-3.3/GenomeAnalysisTK.jar -T IndelRealigner -R ../Reference.fasta -I ../RC3.SAMTOOLSVIEW.bam -targetIntervals ../RC3.GATKREALIGNERTARGETCREATOR.intervals -o ../RC3.GATKINDELREALIGNER.bam 
    

Rk : Change --fix_misencoded_quality_scores -fixMisencodedQuals for --read_filter NotPrimaryAlignment => error :1. Error caching SAM record, ..., which is usually caused by malformed SAM/BAMfiles in which multiple identical copies of a read are present.

samtools index

12- picardtools markDuplicates

  • usr/local/java/latest/bin/java -Xmx12g  -jar /home/sabotf/sources/picard-tools/MarkDuplicates.jar INPUT=../RC3.GATKINDELREALIGNER.bam OUTPUT=../RC3.PICARDTOOLSMARKDUPLICATES.bam METRICS_FILE=../RC3.PICARDTOOLSMARKDUPLICATES.bamDuplicates VALIDATION_STRINGENCY=SILENT

samtools index

_ 13 - GATK

haplotypeCaller

/usr/local/java/latest/bin/java -Xmx12g -jar /home/sabotf/sources/GenomeAnalysisTK/GenomeAnalysisTK.jar -T HaplotypeCaller -R ../Reference.fasta -I ../RC3.PICARDTOOLSMARKDUPLICATES.bam -I ../RC3Single.PICARDTOOLSMARKDUPLICATES.bam -o ../GATKHAPLOTYPECALLER.vcf --filter_mismatching_base_and_quals

/usr/java/jre1.7.0_51/bin/java -Xmx12g -jar /usr/local/GenomeAnalysisTK-3.3/GenomeAnalysisTK.jar -T HaplotypeCaller -R ../Reference.fasta -I ../RC3.PICARDTOOLSMARKDUPLICATES.bam -I ../RC3Single.PICARDTOOLSMARKDUPLICATES.bam -o ../GATKHAPLOTYPECALLERbis.vcf

??Rk: --filter_mismatching_base_and_quals => error= is malformed: BAM file has a read with mismatching number of bases and base qualities.
http://gatkforums.broadinstitute.org/discussion/1429/error-bam-file-has-a-read-with-mismatching-number-of-bases-and-base-qualities

gatk variantFiltration

/usr/local/java/latest/bin/java -Xmx12g -jar /home/sabotf/sources/GenomeAnalysisTK/GenomeAnalysisTK.jar -T VariantFiltration -R ../Reference.fasta --variant ../GATKHAPLOTYPECALLER.vcf -o ../GATKVARIANTFILTRATION.vcf

/usr/java/jre1.7.0_51/bin/java -Xmx12g -jar /usr/local/GenomeAnalysisTK-3.3/GenomeAnalysisTK.jar -T VariantFiltration -R ../Reference.fasta --variant ../GATKHAPLOTYPECALLER.vcf -o ../GATKVARIANTFILTRATION.vcf

gatk selectVariants

/usr/local/java/latest/bin/java -Xmx12g -jar /home/sabotf/sources/GenomeAnalysisTK/GenomeAnalysisTK.jar -T SelectVariants -R ../Reference.fasta --variant ../GATKVARIANTFILTRATION.vcf -o ../GATKSELECTVARIANTS.vcf

/usr/java/jre1.7.0_51/bin/java -Xmx12g -jar /usr/local/GenomeAnalysisTK-3.3/GenomeAnalysisTK.jar -T SelectVariants -R ../Reference.fasta --variant ../GATKVARIANTFILTRATION.vcf -o ../GATKSELECTVARIANTS.vcf


  • samtools view -H bamFile > header.txt 
    printf "@RG\tID:--RG--\tSM:--RG--\tPL:Illumina" >> header.txt
     samtools reheader header.txt $bamFile > reheadedFile