贡献一个SNP/Indel calling pipeline# Biology - 生物学f*h2012-09-20 07:091 楼想在以前的旧ipod上装个skype,但是现在itune里的skype必须要ios7.0以上没法用。请问怎么能装一个老版本的skype?谢谢!
j*p2012-09-20 07:092 楼攒人品,顺便回答一下 iiiir 的问题。我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的predictionpower。从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。以下是GATK 的pipeline假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C_R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中的SNPs (相对于C)假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19Check this website if you have any questions:http://seqanswers.com/wiki/How-to/exome_analysis#step 0 align readsbwa sampe -P -r "@RG\tID:sample\tLB:chr_WGS\tSM:sample\tPL:ILLUMINA" /bwa/indexes/hg19 A_R1.sai A_R2.sai A_R1.fastq A_R2.fastq > A.sambwa sampe -P -r "@RG\tID:control\tLB:chr_WGS\tSM:control\tPL:ILLUMINA" /seq/bwa/indexes/hg19 C_R1.sai C_R2.sai C_R1.fastq C_R2.fastq > C.sam# grep all aligned readsgrep chr A.sam > temp.sammv temp.sam A.samgrep chr C.sam > temp2.sammv temp2.sam C.sam#step 1 run picardjava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinateINPUT=A.sam OUTPUT=A.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=truejava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinateINPUT=C.sam OUTPUT=C.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true#step 2 mark duplicate readsjava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/MarkDuplicates.jar INPUT=A.bam OUTPUT=A.marked.bam METRICS_FILE=metrics CREATE_INDEX=true VALIDATION_STRINGENCY=LENIENTjava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/MarkDuplicates.jar INPUT=C.bam OUTPUT=C.marked.bam METRICS_FILE=metrics CREATE_INDEX=true VALIDATION_STRINGENCY=LENIENT#step 3 run gatk, update gatk installation path if needed, make sure “-Lexonscaptured.bed” is only needed if you are doing ExomeSeq. For WGS, notneeded.java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -R /bwa/indexes/hg19.fa -o A.bam.list -I A.marked.bam -L exonscaptured.bedjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -R /bwa/indexes/hg19.fa -o C.bam.list -I C.marked.bam -L exonscaptured.bed#step 4java -Xmx8g -Djava.io.tmpdir=/tmp -jar /gatk/dist/GenomeAnalysisTK.jar -I A.marked.bam -R /bwa/indexes/hg19.fa -T IndelRealigner -targetIntervals A.bam.list -o A.marked.realigned.bam -L exonscaptured.bedjava -Xmx8g -Djava.io.tmpdir=/tmp -jar /gatk/dist/GenomeAnalysisTK.jar -I C.marked.bam -R /bwa/indexes/hg19.fa -T IndelRealigner -targetIntervals C.bam.list -o C.marked.realigned.bam -L exonscaptured.bed#step 5java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=A.marked.realigned.bam OUTPUT=A.marked.realigned.fixed.bam SO=coordinateVALIDATION_STRINGENCY=LENIENT CREATE_INDEX=truejava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=C.marked.realigned.bam OUTPUT=C.marked.realigned.fixed.bam SO=coordinateVALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true#step 6, dbsnp135.vcf can be downloaded from somewhere online.java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/hg19.fa -knownSites /DBSNP/dbsnp135.vcf -I A.marked.realigned.fixed.bam -TCountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -covCycleCovariate -cov DinucCovariate -recalFile A.recal_data.csv -Lexonscaptured.bedjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/hg19.fa -knownSites /DBSNP/dbsnp135.vcf -I C.marked.realigned.fixed.bam -TCountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -covCycleCovariate -cov DinucCovariate -recalFile C.recal_data.csv -Lexonscaptured.bed#step 7java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/hg19.fa -I A.marked.realigned.fixed.bam -T TableRecalibration --out A.marked.realigned.fixed.recal.bam -recalFile ss1.recal_data.csv -L exonscaptured.bedjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/hg19.fa -I C.marked.realigned.fixed.bam -T TableRecalibration --out C.marked.realigned.fixed.recal.bam -recalFile ss1.recal_data.csv -L exonscaptured.bed#step 8java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -glm BOTH -R /bwa/indexes/hg19.fa -T UnifiedGenotyper -I A.marked.realigned.fixed.recal.bam -I C.marked.realigned.fixed.recal.bam -D /DBSNP/dbsnp135.vcf -o A.snps.AB.vcf -metrics A.snps.metrics -stand_call_conf 50.0 -stand_emit_conf 10.0 -Lexonscaptured.bed#select step, select SNPjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TSelectVariants --variant A.snps.AB.vcf -o A.snpsonly.vcf -selectType SNP -Lexonscaptured.bed#select step2, select Indeljava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TSelectVariants --variant A.snps.AB.vcf -o A.indels.vcf -selectType INDEL -Lexonscaptured.bed#step 9 more filter from HapMap and 1000 Genome projectjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TVariantRecalibrator -input A.snpsonly.vcf -resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.b37.sites.vcf -resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.b37.sites.vcf-an QD -an HaplotypeScore -an MQRankSum -an ReadPosRankSum -an MQ -recalFile A.snps.recal.vcf -tranchesFile A.snp.tranches.vcf -mode SNP -Lexonscaptured.bed#step10java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TApplyRecalibration -input A.snpsonly.vcf -tranchesFile A.snp.tranches.vcf -recalFile A.snps.recal.vcf -o A.snps.filtered.vcf -L exonscaptured.bed#step 11java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TVariantFiltration --variant A.indels.vcf -o A.indels.filtered.vcf --filterExpression "QD < 2.0" --filterExpression "ReadPosRankSum < -20.0" --filterExpression "FS > 200.0" --filterName QDFilter --filterNameReadPosFilter --filterName FSFilter -L exonscaptured.bed#step 12java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TVariantFiltration --variant A.snps.filtered.vcf --clusterSize 3 --clusterWindowSize 10 -o A.snps.finalfilteredAB.vcf -L exonscaptured.bed
a*m2012-09-20 07:094 楼现在call snp已经很成熟了吧。光看你列这些软件就知道了~,C【在 j*p 的大作中提到】: 攒人品,顺便回答一下 iiiir 的问题。: 我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,: SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的prediction: power。: 从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。: 以下是GATK 的pipeline: 假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C: _R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中: 的SNPs (相对于C): 假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19
u*12012-09-20 07:095 楼很感激。我赶紧去对比下你的GATK pipeline和我自己的但问题是,貌似GATK的website更新了,貌似你的pipeline也是V3的吧。。现在都V4了,C【在 j*p 的大作中提到】: 攒人品,顺便回答一下 iiiir 的问题。: 我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,: SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的prediction: power。: 从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。: 以下是GATK 的pipeline: 假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C: _R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中: 的SNPs (相对于C): 假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19
i*r2012-09-20 07:097 楼谢谢,祝你一切顺利!我这种情况怎么搞(目前还是没有弄明白)有2个原始文件,一个是发飙了的SNPs数据,大概70个不同human poupulations,是一个bed文件,结构是这样(右边省略了若干column)chr1 41217 41218 snp 2 + T A dbsnp.108:rs3863625 NNchr1 41255 41256 snp 3 + C T dbsnp.111:rs4543737 NNchr1 41980 41981 snp 4 + A G dbsnp.86:rs806721姑且叫A.bed另一个是个bam文件,别的地方下载的,是另一个human population数据,转成sam文件之后(为了好看把SEQ和QUAL拿掉了):all-hg18_1 0 chr1 39 105 97M * 0 0 () () AS:i:0all-hg18_2 0 chr1 138 114 9M1D86M1D22M * 0 0 () () AS:i:0all-hg18_3 0 chr1 501 254 47M * 0 0 () () AS:i:0姑且叫B.bam我们想知道上面文件里面的SNPs,是否在这个human population也可能存在,所以要对第一个文件(A.bed)里面的每一个SNP位点,看它是不是在第二个文件的区间里面(B.bam),如果是,根据序列和CIGAR找出对应点得碱基。另一个问题是samtools sort后并没有按chromosome排序?$ samtools sort B.bam B.sorted$ samtools view B.sorted.bam | awk '{print $3}' | uniqchr1chr10chr11...chr2chr20...不是我想要的:chr1chr2chr3...
u*d2012-09-20 07:098 楼mark,C【在 j*p 的大作中提到】: 攒人品,顺便回答一下 iiiir 的问题。: 我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,: SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的prediction: power。: 从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。: 以下是GATK 的pipeline: 假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C: _R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中: 的SNPs (相对于C): 假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19
j*p2012-09-20 07:099 楼攒人品,顺便回答一下 iiiir 的问题。我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的predictionpower。从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。以下是GATK 的pipeline假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C_R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中的SNPs (相对于C)假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19Check this website if you have any questions:http://seqanswers.com/wiki/How-to/exome_analysis#step 0 align readsbwa sampe -P -r "@RG\tID:sample\tLB:chr_WGS\tSM:sample\tPL:ILLUMINA" /bwa/indexes/hg19 A_R1.sai A_R2.sai A_R1.fastq A_R2.fastq > A.sambwa sampe -P -r "@RG\tID:control\tLB:chr_WGS\tSM:control\tPL:ILLUMINA" /seq/bwa/indexes/hg19 C_R1.sai C_R2.sai C_R1.fastq C_R2.fastq > C.sam# grep all aligned readsgrep chr A.sam > temp.sammv temp.sam A.samgrep chr C.sam > temp2.sammv temp2.sam C.sam#step 1 run picardjava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinateINPUT=A.sam OUTPUT=A.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=truejava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinateINPUT=C.sam OUTPUT=C.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true#step 2 mark duplicate readsjava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/MarkDuplicates.jar INPUT=A.bam OUTPUT=A.marked.bam METRICS_FILE=metrics CREATE_INDEX=true VALIDATION_STRINGENCY=LENIENTjava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/MarkDuplicates.jar INPUT=C.bam OUTPUT=C.marked.bam METRICS_FILE=metrics CREATE_INDEX=true VALIDATION_STRINGENCY=LENIENT#step 3 run gatk, update gatk installation path if needed, make sure “-Lexonscaptured.bed” is only needed if you are doing ExomeSeq. For WGS, notneeded.java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -R /bwa/indexes/hg19.fa -o A.bam.list -I A.marked.bam -L exonscaptured.bedjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -R /bwa/indexes/hg19.fa -o C.bam.list -I C.marked.bam -L exonscaptured.bed#step 4java -Xmx8g -Djava.io.tmpdir=/tmp -jar /gatk/dist/GenomeAnalysisTK.jar -I A.marked.bam -R /bwa/indexes/hg19.fa -T IndelRealigner -targetIntervals A.bam.list -o A.marked.realigned.bam -L exonscaptured.bedjava -Xmx8g -Djava.io.tmpdir=/tmp -jar /gatk/dist/GenomeAnalysisTK.jar -I C.marked.bam -R /bwa/indexes/hg19.fa -T IndelRealigner -targetIntervals C.bam.list -o C.marked.realigned.bam -L exonscaptured.bed#step 5java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=A.marked.realigned.bam OUTPUT=A.marked.realigned.fixed.bam SO=coordinateVALIDATION_STRINGENCY=LENIENT CREATE_INDEX=truejava -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=C.marked.realigned.bam OUTPUT=C.marked.realigned.fixed.bam SO=coordinateVALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true#step 6, dbsnp135.vcf can be downloaded from somewhere online.java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/hg19.fa -knownSites /DBSNP/dbsnp135.vcf -I A.marked.realigned.fixed.bam -TCountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -covCycleCovariate -cov DinucCovariate -recalFile A.recal_data.csv -Lexonscaptured.bedjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/hg19.fa -knownSites /DBSNP/dbsnp135.vcf -I C.marked.realigned.fixed.bam -TCountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -covCycleCovariate -cov DinucCovariate -recalFile C.recal_data.csv -Lexonscaptured.bed#step 7java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/hg19.fa -I A.marked.realigned.fixed.bam -T TableRecalibration --out A.marked.realigned.fixed.recal.bam -recalFile ss1.recal_data.csv -L exonscaptured.bedjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/hg19.fa -I C.marked.realigned.fixed.bam -T TableRecalibration --out C.marked.realigned.fixed.recal.bam -recalFile ss1.recal_data.csv -L exonscaptured.bed#step 8java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -glm BOTH -R /bwa/indexes/hg19.fa -T UnifiedGenotyper -I A.marked.realigned.fixed.recal.bam -I C.marked.realigned.fixed.recal.bam -D /DBSNP/dbsnp135.vcf -o A.snps.AB.vcf -metrics A.snps.metrics -stand_call_conf 50.0 -stand_emit_conf 10.0 -Lexonscaptured.bed#select step, select SNPjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TSelectVariants --variant A.snps.AB.vcf -o A.snpsonly.vcf -selectType SNP -Lexonscaptured.bed#select step2, select Indeljava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TSelectVariants --variant A.snps.AB.vcf -o A.indels.vcf -selectType INDEL -Lexonscaptured.bed#step 9 more filter from HapMap and 1000 Genome projectjava -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TVariantRecalibrator -input A.snpsonly.vcf -resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.b37.sites.vcf -resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.b37.sites.vcf-an QD -an HaplotypeScore -an MQRankSum -an ReadPosRankSum -an MQ -recalFile A.snps.recal.vcf -tranchesFile A.snp.tranches.vcf -mode SNP -Lexonscaptured.bed#step10java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TApplyRecalibration -input A.snpsonly.vcf -tranchesFile A.snp.tranches.vcf -recalFile A.snps.recal.vcf -o A.snps.filtered.vcf -L exonscaptured.bed#step 11java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TVariantFiltration --variant A.indels.vcf -o A.indels.filtered.vcf --filterExpression "QD < 2.0" --filterExpression "ReadPosRankSum < -20.0" --filterExpression "FS > 200.0" --filterName QDFilter --filterNameReadPosFilter --filterName FSFilter -L exonscaptured.bed#step 12java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -TVariantFiltration --variant A.snps.filtered.vcf --clusterSize 3 --clusterWindowSize 10 -o A.snps.finalfilteredAB.vcf -L exonscaptured.bed
a*m2012-09-20 07:0910 楼现在call snp已经很成熟了吧。光看你列这些软件就知道了~,C【在 j*p 的大作中提到】: 攒人品,顺便回答一下 iiiir 的问题。: 我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,: SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的prediction: power。: 从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。: 以下是GATK 的pipeline: 假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C: _R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中: 的SNPs (相对于C): 假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19
u*12012-09-20 07:0911 楼很感激。我赶紧去对比下你的GATK pipeline和我自己的但问题是,貌似GATK的website更新了,貌似你的pipeline也是V3的吧。。现在都V4了,C【在 j*p 的大作中提到】: 攒人品,顺便回答一下 iiiir 的问题。: 我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,: SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的prediction: power。: 从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。: 以下是GATK 的pipeline: 假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C: _R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中: 的SNPs (相对于C): 假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19
i*r2012-09-20 07:0913 楼谢谢,祝你一切顺利!我这种情况怎么搞(目前还是没有弄明白)有2个原始文件,一个是发飙了的SNPs数据,大概70个不同human poupulations,是一个bed文件,结构是这样(右边省略了若干column)chr1 41217 41218 snp 2 + T A dbsnp.108:rs3863625 NNchr1 41255 41256 snp 3 + C T dbsnp.111:rs4543737 NNchr1 41980 41981 snp 4 + A G dbsnp.86:rs806721姑且叫A.bed另一个是个bam文件,别的地方下载的,是另一个human population数据,转成sam文件之后(为了好看把SEQ和QUAL拿掉了):all-hg18_1 0 chr1 39 105 97M * 0 0 () () AS:i:0all-hg18_2 0 chr1 138 114 9M1D86M1D22M * 0 0 () () AS:i:0all-hg18_3 0 chr1 501 254 47M * 0 0 () () AS:i:0姑且叫B.bam我们想知道上面文件里面的SNPs,是否在这个human population也可能存在,所以要对第一个文件(A.bed)里面的每一个SNP位点,看它是不是在第二个文件的区间里面(B.bam),如果是,根据序列和CIGAR找出对应点得碱基。另一个问题是samtools sort后并没有按chromosome排序?$ samtools sort B.bam B.sorted$ samtools view B.sorted.bam | awk '{print $3}' | uniqchr1chr10chr11...chr2chr20...不是我想要的:chr1chr2chr3...
u*d2012-09-20 07:0914 楼mark,C【在 j*p 的大作中提到】: 攒人品,顺便回答一下 iiiir 的问题。: 我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,: SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的prediction: power。: 从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。: 以下是GATK 的pipeline: 假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C: _R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中: 的SNPs (相对于C): 假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19
s*r2012-09-20 07:0915 楼请问有lz用过GATK最新的HaplotypeCaller么GATK自己说这是他们最先进的variants caller不过我试了下 发现很多genotyps都和我在IGV里看的矛盾,C【在 j*p 的大作中提到】: 攒人品,顺便回答一下 iiiir 的问题。: 我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,: SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的prediction: power。: 从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。: 以下是GATK 的pipeline: 假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C: _R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中: 的SNPs (相对于C): 假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19
u*12012-09-20 07:0917 楼你要确定你BWA产生bam file时候使用的ref,和你GATK command里的ref,是同一个reference;GATK对这个要求很严,尤其不同版本的ref,比如g1k_37和hg19;g1k会包含一些无法assembly的contig。总之ref里的chromosome的多少,排序都很重要,同时要注意prefix比如“chr”;20和chr20是不一样的,GATK没办法识别的【在 x*****d 的大作中提到】: 我用的是Tophat1.4加freebayes。 GATK老是说我的reference有问题。
x*d2012-09-20 07:0918 楼多谢指点!我用的都是Tophat主页上面的hg19和mm9 reference。下次再试试。【在 u*********1 的大作中提到】: 你要确定你BWA产生bam file时候使用的ref,和你GATK command里的ref,是同一个: reference;GATK对这个要求很严,尤其不同版本的ref,比如g1k_37和hg19;g1k会包: 含一些无法assembly的contig。总之ref里的chromosome的多少,排序都很重要,同时: 要注意prefix比如“chr”;20和chr20是不一样的,GATK没办法识别的
x*d2012-09-20 07:0919 楼楼主少了一步吧?bwa aln -t 4 -f -I fastq>bwa aln -t 4 -f -I fastq>,C【在 j*p 的大作中提到】: 攒人品,顺便回答一下 iiiir 的问题。: 我们尝试过好几种不同的SNP calling的方法,包括GATK, Samtools, Varscan,: SeqGenes, 等,并且做了SNP array 作为gold standard比较各种方法的prediction: power。: 从我们的经验,BWA + GATK 最好,sensitivity 和 specificity 都在95%以上。: 以下是GATK 的pipeline: 假设你有一个control 样品C 和一个样本样品A的pair-end sequencing,共4个文件,C: _R1.fastq, C_R2.fastq, A_R1.fastq and A_R2.fastq如何通过BWA/GATK去找样品A中: 的SNPs (相对于C): 假设assembly 用的是hg19,你的BWA index 在这里:/bwa/indexes/hg19
g*g2012-09-20 07:0920 楼For SNP calling, there is a better one came out recently:http://www.nature.com/ncomms/journal/v3/n12/abs/ncomms2256.html
x*d2012-09-20 07:0921 楼看了,这个流程要求pileup.问题是我的样品是100 million PE read. pileup后文件实在太大了,没有GATK占的空间少.【在 g*****g 的大作中提到】: For SNP calling, there is a better one came out recently:: http://www.nature.com/ncomms/journal/v3/n12/abs/ncomms2256.html