大家可以根据我写的千人基因组计划ftp资源介绍,找到下面3个样本的测序文件:
然后走一遍外显子数据分析流程!
一个简单的流程脚本如下:
[Shell] 纯文本查看 复制代码 ################################
#This pipeline is for single sample variants detection including
#Single Nucleotide Variants and Short Indels from
#Whole Genome Sequencing (WGS) or Whole Exome Sequencing (WES)
#with a minimum 10x of average depth of coverage.
#
#Sequencing platform: illumina HiSeq/MiSeq
#Pair-end reads with length from 50bp to 150bp.
################################
#env settings, Java Runtime required. (openjdk-7 is used here)
CUROVERSE_HOME=/mnt/curoverse
APP=$CUROVERSE_HOME/app
DATA=$CUROVERSE_HOME/data
################################
#prepare applications
################################
#1. BWA
wget [url]http://downloads.sourceforge.net/project/bio-bwa/bwa-0.7.9a.tar.bz2[/url]
tar -jxvf bwa-0.7.9a.tar.bz2
cd bwa-0.7.9a && make && cd $CUROVERSE_HOME
mkdir -p $APP/bwa/0.7.9a/
find bwa-0.7.9a -executable -type f -print0 | xargs -0 -I {} cp {} $APP/bwa/0.7.9a/
rm -fr bwa*
#2. samtools
wget [url]http://downloads.sourceforge.net/project/samtools/samtools/0.1.19/samtools-0.1.19.tar.bz2[/url]
tar -jxvf samtools-0.1.19.tar.bz2
cd samtools-0.1.19 && make && cd $CUROVERSE_HOME
mkdir -p $APP/samtools/0.1.19/
find samtools-0.1.19 -executable -type f -print0 | xargs -0 -I {} mv {} $APP/samtools/0.1.19/
rm -fr samtools*
#3. picard
wget [url]http://downloads.sourceforge.net/project/picard/picard-tools/1.114/picard-tools-1.114.zip[/url]
unzip picard-tools-1.114.zip
mkdir -p $APP/picard/1.114/
mv picard-tools-1.114/* $APP/picard/1.114/
rm -fr picard*
#4. GATK
#no automatic way to install GATK.
Register and download the GATK package. unpack it to /mnt/curoverse/app/gatk/3.1-1/
################################
#prepare reference data
################################
mkdir -p $DATA/fasta/hg19/ && cd $DATA/fasta/hg19/
for i in $(seq 1 22) X Y M; do echo $i; wget [url]http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr[/url]${i}.fa.gz; done
gunzip *.gz
for i in $(seq 1 22) X Y M; do cat chr${i}.fa >> hg19.fasta; done
rm -fr chr*.fasta
#create index
$APP/samtools/0.1.19/samtools faidx hg19.fasta
#create dictionary
java -Xmx8g -jar $APP/picard/0.1.19/CreateSequenceDictionary.jar R=hg19.fasta O=hg19.dict
java -Xmx8g -jar /home/hadoop/curoverse/app/picard/1.114/CreateSequenceDictionary.jar R=hg19.fasta O=hg19.dict
cd $CUROVERSE_HOME
################################
#prepare library files (dbSNP, 1000G, etc)
################################
mkdir -p $DATA/lib/hg19/ && cd $DATA/lib/hg19/
#dbSNP138
#wget [url]http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp138Common.txt.gz[/url]
#GATK resource bundle
GATK_FTP=ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/2.8/hg19
#download library file for GATK's VariantRecalibrator
for f in hapmap_3.3.hg19.vcf.gz 1000G_omni2.5.hg19.vcf.gz 1000G_phase1.indels.hg19.vcf.gz Mills_and_1000G_gold_standard.indels.hg19.vcf.gz
do
curl --remote-name ${GATK_FTP}/${f}
done
################################
#Demo sample data (NA12878 from 1000G)
#NA12878: pair-end, WGS, 50x depth, Illumina HiSeq 2000 system
#URL: [url]http://www.ebi.ac.uk/ena/data/view/ERX237515[/url]
#NA12878_R1: [url]ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_1.fastq.gz[/url]
#NA12878_R2: [url]ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_2.fastq.gz[/url]
################################
#~100GB disk space required to download the NA12878
wget [url]ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_1.fastq.gz[/url]
wget [url]ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_2.fastq.gz[/url]
################################
#Begin of pipeline
#For testing purpose, two PE fastq files R1 and R2 were generated by pIRS ([url]http://code.google.com/p/pirs/[/url])
#and only chr1 of hg19 was used as reference genome
cd /mnt/
#0.1 create simulated PE Reads on 50x coverage
pirs simulate -i $DATA/hg19.chr1/hg19.chr1.fa \
-s ~/app/pIRS_111/Profiles/Base-Calling_Profiles/humNew.PE100.matrix.gz \
-d ~/app/pIRS_111/Profiles/GC-depth_Profiles/humNew.gcdep_200.dat \
-b ~/app/pIRS_111/Profiles/InDel_Profiles/phixv2.InDel.matrix \
-l 100 -x 50 -Q 33 -o A
R1=/mnt/A_100_500_1.fq.gz
R2=/mnt/A_100_500_2.fq.gz
#0.2 build genome index for BWA
$APP/app/bwa/0.7.9a/bwa index -p hg19.chr1 chr1.fa
#1. align with BWA, with 4 threads (-t 4). Reading Group is required (-R ...).
time $APP/bwa/0.7.9a/bwa mem -t 4 -R '@RG\tID:group_id\tPL:illumina\tSM:sample_id' $DATA/hg19.chr1/hg19.chr1 ${R1} ${R2} > A.chr1.sam
#2. Sort into BAM
java -Xmx4g -Djava.io.tmpdir=/tmp -jar $APP/picard/1.114/SortSam.jar \
CREATE_INDEX=True SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT \
INPUT=A.chr1.sam OUTPUT=A.chr1.sort.bam
#3. remove PCR duplicate
java -Xmx4g -Djava.io.tmpdir=/tmp -jar $APP/picard/1.114/MarkDuplicates.jar \
CREATE_INDEX=true REMOVE_DUPLICATES=True ASSUME_SORTED=True VALIDATION_STRINGENCY=LENIENT METRICS_FILE=/dev/null \
INPUT=A.chr1.sort.bam OUTPUT=A.chr1.sort.dedup.bam
#4. fix mate information
java -Djava.io.tmpdir=/tmp/ -jar $APP/picard/1.114/FixMateInformation.jar \
SO=coordinate VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true \
INPUT=A.chr1.sort.dedup.bam OUTPUT=A.chr1.sort.dedup.mate.bam
#5. local alignment around Indels. create de novo intervals
java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T RealignerTargetCreator \
-I A.chr1.sort.dedup.mate.bam -o A.intervals
java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T IndelRealigner \
-targetIntervals A.intervals \
-I A.chr1.sort.dedup.mate.bam -o A.chr1.sort.dedup.mate.relaign.bam
#6. base recalibrator, use dbSNP139 as "truth". If we know the population of subject sample, we can
#use specific SNP datasets extract from 1000G with same population.
java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T BaseRecalibrator \
-knownSites $DATA/dbsnp_138.hg19.vcf \
-o A.recal \
-I A.chr1.sort.dedup.mate.relaign.bam
java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T PrintReads \
-BQSR A.recal \
-o A.chr1.sort.dedup.mate.relaign.recal.bam \
-I A.chr1.sort.dedup.mate.relaign.bam
#7. call variants with GATK's UnifiedGenotyper (faster, lots of raw calls, must run filter after)
time java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T UnifiedGenotyper \
-glm BOTH \
-D $DATA/dbsnp_138.hg19.vcf \
-metrics A.snps.metrics \
-stand_call_conf 50.0 -stand_emit_conf 10.0 \
-o A.ug.vcf \
-I A.chr1.sort.dedup.mate.relaign.recal.bam
#8. call variants with GATK's HaplotypeCaller (2x slower than UG, much less accurate calls)
time java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T HaplotypeCaller \
--dbsnp $DATA/dbsnp_138.hg19.vcf \
-stand_call_conf 50.0 -stand_emit_conf 10.0 \
-o A.hc.vcf \
-I A.chr1.sort.dedup.mate.relaign.recal.bam
#End of pipeline
################################
|