Alice:Exome Sequencing Analysis/analysis pipeline
Jump to navigation
Jump to search
Update[edit]
- Since Athurva updated the GATK to v3.3 back in 2015 (also replaced UnifiedGenotyper with HaplotypeCaller), so I am going to use the new version to run analysis on mouse samples
- here are the software location:
$mapper = "/home/kunzhang/softwares/bwa-0.7.5a/bwa"; $samtool = "/home/kunzhang/softwares/samtools-latest/samtools"; $picardDir="java -Xmx4g -jar /home/kunzhang/softwares/picard-tools-latest/"; $gatk_dir="/home/kunzhang/softwares/GenomeAnalysisTK-3.3/";
- These are the location of the old files:
mouse files: $dbsnp="/media/Ext12T/GenomeDB/MmGenome/mm9/GATK_1.2/mm9_dbSNP132.vcf"; $CCDS_bed="/media/Ext12T/GenomeDB/MmGenome/mm9/GATK_1.2/CCDS.20120710.UCSC.chr.bed"; $template = "/media/Ext12T/GenomeDB/MmGenome/mm9/GATK_1.2/Mm9.fasta"; $fai_file = "/media/Ext12T/GenomeDB/MmGenome/mm9/GATK_1.2/Mm9.fasta.fai"; human files: $dbsnp="/home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/dbsnp_132.b37.vcf"; $CCDS_bed="/home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/CCDS.20120105.UCSC.nochr.bed"; $template = "/home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/human_g1k_v37.fasta"; $fai_file = "/home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/human_g1k_v37.fasta.fai";
- these are the location of new files:
For human sample: $dbsnp="/home/kunzhang/HsGenome/resources_v2.8_b37/dbsnp_138.b37.vcf"; $CCDS_bed="/home/kunzhang/HsGenome/resources_v2.8_b37/CCDS.20120105.UCSC.nochr.bed"; $known_indels = "/home/kunzhang/HsGenome/resources_v2.8_b37/Mills_and_1000G_gold_standard.indels.b37.vcf"; $template = "/home/kunzhang/HsGenome/resources_v2.8_b37/human_g1k_v37_decoy.fasta"; $fai_file = "/home/kunzhang/HsGenome/resources_v2.8_b37/human_g1k_v37_decoy.fasta.fai"; For mouse sample: $dbsnp="/media/Ext12T/GenomeDB/MmGenome/mm9/GATK_1.2/mm9_dbSNP132.vcf"; $CCDS_bed="/media/Ext12T/GenomeDB/MmGenome/mm9/GATK_1.2/CCDS.20120710.UCSC.chr.bed"; $template = "/home/zhl002/genomeDB/mm9/Mm9.fasta"; $fai_file = "/home/zhl002/genomeDB/mm9/Mm9.fasta.fai"
- edit the updated pipeline with the files for mm9 listed above, also there is no mouse vcf file for indels
Human/Mouse Genome Resequencing Pipeline[edit]
- Everything below is adapted and modified from Athurva's wiki page for my reference: Link
- newest version of the pipeline:
- genome-miner:
- Human: /home/kunzhang/bin/variantCallerBwaGATK_05012015.pl
- Mouse: /home/zhl002/bin/variantCallerBwaGATK_mouse_2015.pl
- genome-miner:
- tscc:
- Human: /home/zhl002/bin/variantCallerBwaGATK_05012015.pl
- Mouse: /home/zhl002/bin/variantCallerBwaGATK_mouse_2015.pl (debugged and working)
- tscc:
- last version of the pipeline:
- genome-miner:
- Human: /home/ajgore/AG_Ext12T/GATK_01022012/variantCallerBwaGATK-latest/variantCallerBwaGATK.pl
- Mouse: /home/ajgore/AG_Ext12T/GATK_01022012/variantCallerBwaGATK-latest/variantCallerBwaGATK_mouse.pl
- genome-miner:
Requirements[edit]
- The below paths list genome-miner locations for human files.
- NCBI Human Genome Reference Sequence, v37 (/home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/human_g1k_v37.fasta)
- CCDS BED file (/home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/CCDS.20120105.UCSC.nochr.bed)
- dbSNP in VCF format (/home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/dbsnp_132.b37.vcf)
- GATK from 01022012 (/home/ajgore/AG_Ext12T/GATK_01022012/GenomeAnalysisTK-1.4-1-g1b6d528/)
- BWA 0.5.8 or later (/home/kunzhang/softwares/bwa-latest)
- Samtools 0.1.7 or later (/home/kunzhang/softwares/samtools-latest)
- Picard-tools 1.38 or later (/home/kunzhang/softwares/picard-tools-latest)
Instructions[edit]
- For mouse data, follow procedure as below, but append "_mouse.pl" to variantCallerBwaGATK command.
Generate data annotation file (.info file)[edit]
- First, generate a text file containing your data information. This file should look like the following:
genome-miner 1 exome SampleID LibraryID FlowCell Lane Barcode File1 File2 CV-iF AL-Exome-CV-iF-PP HL022 s5 NA /media/Syn_15T/Ext9T_Backup/SeqStore2009/091109_HL022/s_5_sequence.txt.gz CV-iF AL-Exome-CV-iF-PP HL022 s6 NA /media/Syn_15T/Ext9T_Backup/SeqStore2009/091109_HL022/s_6_sequence.txt.gz CV-iF AL-Exome-CV-iF-PP HL022 s7 NA /media/Syn_15T/Ext9T_Backup/SeqStore2009/091109_HL022/s_7_sequence.txt.gz CV-iF AL-Exome-CV-iF-PP HL025 s7 NA /media/Syn_15T/Ext9T_Backup/SeqStore2009/091123_HL025/s_7_sequence.txt.gz CV-iF AL-Exome-CV-iF-PP HL026 s5 NA /media/Syn_15T/Ext9T_Backup/SeqStore2009/091127_HL026/s_5_sequence.txt.gz CV-iF AL-Exome-CV-iF-NG HL053 s2 NA /media/SeqStore2/100503_HL053/s_2_1_sequence.txt.gz /media/SeqStore2/100503_HL053/s_2_2_sequence.txt.gz
- This file is a tab-delimited text file containing 7 columns
- The first three lines list script parameters.
- Line 1: Please write either "genome-miner" or "tscc" based on the system of choice
- Line 2: Please write "1" to remove clonal reads or "0" to keep clonal reads
- Line 3: Please write "exome" to limit realignment and variant calling to coding regions, or "genome" to look at all regions
- Line 4: MUST contain the header line listed above.
- Each line represents one SE fastq file or a pair of PE fastq files.
- Detailed description of each column:
- Column 1: SampleID
- This column should contain a short identifier describing your sample(s).
- If multiple samples are provided, a multi-sample BAM file will be generated and multi-sample variant calling will be performed.
- Column 2: LibraryID
- This column should describe details of your library.
- Generally, you should include the initials of the person making the library, the type (Exome, Genome, BSPP, etc), and the sampleID.
- If you generated multiple sequencing libraries for your samples, please provide them unique identifiers here to facilitate proper removal of clonal reads
- For example, here I have denoted that one library was made using padlock probes and one was made using NimbleGen Exome Capture.
- Column 3: FlowCell
- This column should contain the flow cell in which the fastq file was sequenced.
- Column 4: Lane
- This column should contain the letter "s" followed by the lane number.
- Column 5: Barcode
- This column should list the barcode used for the library.
- Examples: "NA" for no barcode or "Indx7" for barcode #7, etc.
- Column 6: First End FASTQ file
- This file should list the full path to the first end FASTQ file.
- For Single-End data, simply provide the path to the FASTQ file.
- Column 7: Second End FASTQ File
- This column should list the full path to the second end FASTQ file
- For single-end data, this column should be left blank.
- Column 1: SampleID
- Save this file with a .info extension. I will call the file "prefix.info" in this example.
- Upload this file to genome-miner to your folder of choice. For this example, we will use /home/ajgore/TestRun.
Run pipeline (Genome-Miner)[edit]
- To run the pipeline, issue the following terminal command from the folder you uploaded the .info file to (in this case /home/ajgore/TestRun):
nohup /home/ajgore/AG_Ext12T/GATK_01022012/variantCallerBwaGATK-latest/variantCallerBwaGATK.pl prefix.info &> prefix.status &
- You can monitor the status of the pipeline by looking at the "prefix.status" file using less:
less prefix.status
- The output of the pipeline will be two files:
- prefix.fixed.bam
- prefix.snp.raw.vcf
Run pipeline (Triton)[edit]
- On Triton, you will need to make a job file that submits the pipeline to the cluster.
- A template job file can be found at: triton-login:/projects/zhang-lab/softwares/GenomeAnalysisTK_01022012/variantCallerBwaGATK_03282012/prefix.job
- For this case, the job file (named prefix.job) should look like:
#!/bin/csh #PBS -q small #PBS -l walltime=24:00:00 #PBS -o variantCallerBwaGATK.outputfile.txt #PBS -e variantCallerBwaGATK.errorfile.txt #PBS -N prefix.info #PBS -V #PBS -M ajgore@eng.ucsd.edu #PBS -m abe #PBS -A zhang-lab cd /home/ajgore/TestRun /projects/zhang-lab/softwares/GenomeAnalysisTK_01022012/variantCallerBwaGATK-latest/variantCallerBwaGATK.pl prefix.info >& prefix.status
- NOTE: The ampersand in the pipeline command has been moved after the redirect due to differences between bash and csh.
- To run the pipeline, issue the following terminal command from the folder you uploaded the .info file to (in this case /home/ajgore/TestRun):
qsub prefix.job
- You can monitor the status of the pipeline by looking at the "prefix.status" file using less:
less prefix.status
- The output of the pipeline will be two files:
- prefix.fixed.bam
- prefix.snp.raw.vcf
Exome Variant Calling Pipeline[edit]
- Pipeline can be found in:
- genome-miner:
- Human: /home/ajgore/AG_Ext12T/GATK_01022012/variantFiltration-latest/variantFiltration.pl
- Mouse: /home/ajgore/AG_Ext12T/GATK_01022012/variantFiltration-latest/variantFiltration_mouse.pl
- triton:
- Human: /projects/zhang-lab/softwares/GenomeAnalysisTK_01022012/variantFiltration-latest/variantFiltration.pl
- Mouse: /projects/zhang-lab/softwares/GenomeAnalysisTK_01022012/variantFiltration-latest/variantFiltration_mouse.pl
- This link will always point to the updated version
- genome-miner:
- Updated on 07/11/2012
Requirements[edit]
- All required files are currently on genome-miner and triton cluster. The below paths list genome-miner locations.
- GATK from 01022012 (/home/ajgore/AG_Ext12T/GATK_01022012/GenomeAnalysisTK-1.4-1-g1b6d528/)
- Samtools 0.1.7 or later (/home/kunzhang/softwares/samtools-latest)
- NCBI Human Genome Reference Sequence, v37 (/home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/human_g1k_v37.fasta)
- Note: Due to bug in samtools pileup, needed to create a separate index file on a symlinked genome.
- This can be found here (also present on triton): /home/ajgore/AG_Ext12T/GATK_01022012/resources_1.2/human_g1k_v37_samtools7.fasta
Instructions[edit]
Gather necessary files[edit]
- Two output files from mapping pipeline are required:
- .fixed.bam file
- .snp.raw.vcf file
- Either run the pipeline in the same directory as these files or use the full paths to these files in the command below.
Variant filtration pipeline[edit]
- To run the pipeline, issue the following terminal command from the folder you uploaded the .info file to (in this case /home/ajgore/TestRun):
nohup /home/ajgore/AG_Ext12T/GATK_01022012/variantFiltration-latest/variantFiltration.pl prefix.snp.raw.vcf prefix.fixed.bam genome-miner &> prefix.filter.status &
- The first argument should be the .vcf file, the second argument should be the .fixed.bam file, and the third argument should be the system (either genome-miner or triton)
- You can monitor the status of the pipeline by looking at the "prefix.status" file using less:
less prefix.filter.status
- The output of the pipeline will be five files:
- prefix.ind.filtered.vcf - VCF file containing filtered indel calls both passing and failing filters (GATK Best Practices Exome v3 Filters)
- prefix.snp.filtered.GATK_New.vcf - VCF file containing snp calls both passing and failing filters (GATK Best Practices Exome v3 Filters)
- prefix.snp.GATK_New.pileup - Pileup file containing filtered snp calls ONLY (GATK Best Practices Exome v3 Filters)
- prefix.snp.filtered.GATK_Old.vcf - VCF file containing snp calls both passing and failing filters (GATK Best Practices Exome v1 Filters)
- prefix.snp.GATK_Old.pileup - Pileup file containing filtered snp calls ONLY (GATK Best Practices Exome v1 Filters)
variant candidate filtering[edit]
- location of scripts on genome-miner:
/home/zhl002/bin/getPileupFromVCF.pl $progenitor(/contol).fixed.bam $sample.snp.filtered.fixed.vcf > $control.$sample.pileup /home/zhl002/bin/getPileupFromVCF.pl $sample.fixed.bam $sample.snp.filtered.fixed.vcf > $sample.$sample.pileup /home/zhl002/bin/pileup2variantsPair.pl $control.$sample.pileup $sample.$sample.pileup > $sample.snc grep GOH $sample.snc | grep -v rs | grep "=0" > $sample.snc.filtered
- first two lines make a limited pileup file for both the comtrol and the sample
- then you do the subtraction and apply the filter
- filters listed above as example:
- only look at gains of heterozygosity (grep GOH),
- only look at non-dbsnp sites (grep -v rs)
- only consider cases where we don't see the minor allele in the control
- the other filters (chi-sq test, etc) are done inside the perl script, but leave these three outside so we can look at the unfiltered .snc file and the filtered