ILRI Research Computing

User Tools

Site Tools

You are here: start » mkatari-bioinformatics-august-2013-gatknotes
mkatari-bioinformatics-august-2013-gatknotes

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
Next revisionBoth sides next revision
mkatari-bioinformatics-august-2013-gatknotes [2014/06/11 12:52] – mkatarimkatari-bioinformatics-august-2013-gatknotes [2015/06/08 07:45] – mkatari
Line 29: Line 29:
   * create a bam index for dedup bam files   * create a bam index for dedup bam files
   * realign samples   * realign samples
 +
 +Once they are all processed
   * merge all samples to one bam file   * merge all samples to one bam file
   * sort and index this merged bam file   * sort and index this merged bam file
Line 36: Line 38:
 bowtie2 -x PTC_Human -U Cohen.fastq -S Cohen.sam bowtie2 -x PTC_Human -U Cohen.fastq -S Cohen.sam
 samtools view -bS Cohen.sam > Cohen.bam samtools view -bS Cohen.sam > Cohen.bam
-bowtie2 -x PTC_Human -U Sherman.fastq -S Sherman.sam 
-samtools view -bS Sherman.sam > Sherman.bam 
 </code> </code>
  
-The picard method to sort is preferred by GATK+The picard method to sort is preferred by GATK. In some cases PICARD uses the temp directory to do its sorting. You may run into an error that complains about running out of space. To avoid this problem simply create your own tmp directory and tell java that it should use it. See details [[https://www.biostars.org/p/42613/|here]]. 
 <code> <code>
-java -jar /export/apps/picard-tools/1.112/SortSam.jar \+mkdir /var/scratch/mkatari 
 +mkdir /var/scratch/mkatari/tmp 
 + 
 +java -Djava.io.tmpdir=/var/scratch/mkatari/tmp -jar /export/apps/picard-tools/1.112/SortSam.jar \
    INPUT=Cohen.bam \    INPUT=Cohen.bam \
    OUTPUT=Cohen.sorted.bam \    OUTPUT=Cohen.sorted.bam \
    SORT_ORDER=coordinate    SORT_ORDER=coordinate
        
-java -jar /export/apps/picard-tools/1.112/SortSam.jar \ +
-   INPUT=Sherman.bam \ +
-   OUTPUT=Sherman.sorted.bam \ +
-   SORT_ORDER=coordinate+
 </code> </code>
  
Line 56: Line 57:
  
 <code> <code>
-java -jar /export/apps/picard-tools/1.112/AddOrReplaceReadGroups.jar \ 
-   INPUT=Sherman.sorted.bam \ 
-   OUTPUT=ShermanRG.bam \ 
-   RGLB=Sherman \ 
-   RGPL=IonTorrent \ 
-   RGPU=None \ 
-   RGSM=Sherman 
  
-java -jar /export/apps/picard-tools/1.112/AddOrReplaceReadGroups.jar \+java -Djava.io.tmpdir=/var/scratch/mkatari/tmp -jar /export/apps/picard-tools/1.112/AddOrReplaceReadGroups.jar \
    INPUT=Cohen.sorted.bam \    INPUT=Cohen.sorted.bam \
    OUTPUT=CohenRG.bam \    OUTPUT=CohenRG.bam \
Line 75: Line 69:
 This will remove any reads that map to the same exact place. It is helpful to get rid of artifacts.  This will remove any reads that map to the same exact place. It is helpful to get rid of artifacts. 
 <code> <code>
-java -jar /export/apps/picard-tools/1.112/MarkDuplicates.jar \+ 
 +java -Djava.io.tmpdir=/var/scratch/mkatari/tmp -jar /export/apps/picard-tools/1.112/MarkDuplicates.jar \
    INPUT=CohenRG.bam \    INPUT=CohenRG.bam \
    OUTPUT=Cohen.dedup.bam \    OUTPUT=Cohen.dedup.bam \
    METRICS_FILE=Cohen.dedup.metrics \    METRICS_FILE=Cohen.dedup.metrics \
    REMOVE_DUPLICATES=TRUE \    REMOVE_DUPLICATES=TRUE \
-   ASSUME_SORTED=TRUE+   ASSUME_SORTED=TRUE \ 
 +   MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000
  
-java -jar /export/apps/picard-tools/1.112/MarkDuplicates.jar \ 
-   INPUT=ShermanRG.bam \ 
-   OUTPUT=Sherman.dedup.bam \ 
-   METRICS_FILE=Sherman.dedup.metrics \ 
-   REMOVE_DUPLICATES=TRUE \ 
-   ASSUME_SORTED=TRUE 
 </code> </code>
  
Line 93: Line 83:
 <code> <code>
 samtools index Cohen.dedup.bam  samtools index Cohen.dedup.bam 
-samtools index Sherman.dedup.bam  
  
 #identifying indels #identifying indels
-java -Xmx2g -jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \+java -Xmx2g -Djava.io.tmpdir=/var/scratch/mkatari/tmp -jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \
    -T RealignerTargetCreator \    -T RealignerTargetCreator \
    -R PTC_Human.fasta \    -R PTC_Human.fasta \
Line 102: Line 91:
    -o CohenforIndelRealigner.intervals    -o CohenforIndelRealigner.intervals
    
- #identifying indels 
-java -Xmx2g -jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \ 
-   -T RealignerTargetCreator \ 
-   -R PTC_Human.fasta \ 
-   -I Sherman.dedup.bam \ 
-   -o ShermanforIndelRealigner.intervals 
  
    
- java -Xmx4g -jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \+ java -Xmx4g -Djava.io.tmpdir=/var/scratch/mkatari/tmp -jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \
    -T IndelRealigner \    -T IndelRealigner \
    -R PTC_Human.fasta \    -R PTC_Human.fasta \
Line 116: Line 99:
   -targetIntervals CohenforIndelRealigner.intervals \   -targetIntervals CohenforIndelRealigner.intervals \
    -o Cohen.dedup.realign.bam    -o Cohen.dedup.realign.bam
- 
- java -Xmx4g -jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \ 
-   -T IndelRealigner \ 
-   -R PTC_Human.fasta \ 
-   -I Sherman.dedup.bam \ 
-  -targetIntervals ShermanforIndelRealigner.intervals \ 
-   -o Sherman.dedup.realign.bam 
  
 </code> </code>
  
-Now we merge the bam files and then sort and index them+In some cases there may be a need to clean the sam/bam file(s) (soft-trimming the coordinates). To do this use CleanSam in Picard tools. You may want to just do it to all to avoid the error in a workflow, but it may not be necessary.
  
 <code> <code>
-java -jar /export/apps/picard-tools/1.112/MergeSamFiles.jar \+java -Djava.io.tmpdir=/var/scratch/mkatari/tmp -jar /export/apps/picard-tools/1.112/CleanSam.jar \
    INPUT=Sherman.dedup.realign.bam \    INPUT=Sherman.dedup.realign.bam \
 +   OUTPUT=Sherman.clean.dedup.realign.bam
 +</code>
 +
 +Now we merge the bam files and then sort and index them. If you cleaned the bam file, remember to use the cleaned ones.
 +
 +<code>
 +java -Djava.io.tmpdir=/var/scratch/mkatari/tmp -jar /export/apps/picard-tools/1.112/MergeSamFiles.jar \
 +   INPUT=Sherman.clean.dedup.realign.bam \
    INPUT=Cohen.dedup.realign.bam \    INPUT=Cohen.dedup.realign.bam \
    OUTPUT=ShermanCohenMerged.bam      OUTPUT=ShermanCohenMerged.bam  
Line 140: Line 124:
  
  
-Finall !! run gatk+Finally !! run gatk
  
 <code> <code>
-java -jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \+java -Djava.io.tmpdir=/var/scratch/mkatari/tmp -jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \
    -T UnifiedGenotyper \    -T UnifiedGenotyper \
    -I ShermanCohenMerged.sorted.bam \    -I ShermanCohenMerged.sorted.bam \
Line 151: Line 135:
    -glm SNP \    -glm SNP \
    -o PTC_human.gatk.vcf    -o PTC_human.gatk.vcf
 +
 +</code>
 +
 +If you want to load the vcf file into IGV, remember to index it first.
 +<code>
 +module load igvtools
 +igvtools index PTC_human.gatk.vcf
 +</code>
 +
 +If you would like to generate a table of from the vcf file use the following command
 +<code>
 +java --Djava.io.tmpdir=/var/scratch/mkatari/tmp jar /export/apps/GenomeAnalysisTK/GenomeAnalysisTK-2.3-9-ge5ebf34/GenomeAnalysisTK.jar \
 +     -R PTC_Human.fasta
 +     -T VariantsToTable \
 +     -V PTC_human.gatk.vcf \
 +     -F CHROM -F POS -F ID -F QUAL -F AC \
 +     -GF GT -GF GQ \
 +     -o PTC_human.gatk.vcf.table
 +</code>
 +
 +In order to filter your vcf file based on quality measures, depth, and also statistical significance, you can use variant filter option in the gatk toolkit. Below is an example of suggested filters for data that has low coverage.
 +
 +<code>
 +java -Xmx2g -jar /export/apps/gatk/3.1.1/GenomeAnalysisTK.jar \
 +    -R /home/emasumba/cassavaV5_0.chromsomesRomanNumerals.fa \
 +    -T VariantFiltration \
 +    -o namikonga_albert_filter.vcf \
 +    --variant /home/emasumba/Namikonga_Albert.gatk.vcf \
 +    --filterExpression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || HaplotypeScore >13.0" \
 +    --filterName "mannyfilter"
 +
 +</code>
 +
 +Good descriptions of the different information on vcf files [[https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_annotator_HaplotypeScore.php|GATK Docs]]
 +
 +Finally to save the SNPs that passed your filter, you simply use the selectvariant tool.
 +
 +<code>
 +
 +java -Xmx2g -jar /export/apps/gatk/3.1.1/GenomeAnalysisTK.jar \
 +    -T SelectVariants \
 +    --variant namikonga_albert_filter.vcf \
 +    -o namikonga_albert_filter_only.vcf \
 +    -ef \
 +    -R /home/emasumba/cassavaV5_0.chromsomesRomanNumerals.fa
 +
 </code> </code>
mkatari-bioinformatics-august-2013-gatknotes.txt · Last modified: 2016/08/17 08:37 by mkatari