贡献一个SNP/Indel calling pipeline# Biology - 生物学
f*h
1 楼
想在以前的旧ipod上装个skype，但是现在itune里的skype必须要ios7.0以上没法用。
请问怎么能装一个老版本的skype？谢谢！
请问怎么能装一个老版本的skype？谢谢！
j*p
2 楼
攒人品，顺便回答一下 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
Check this website if you have any questions:
http://seqanswers.com/wiki/How-to/exome_analysis
#step 0 align reads
bwa 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.sam
bwa 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 reads
grep chr A.sam > temp.sam
mv temp.sam A.sam
grep chr C.sam > temp2.sam
mv temp2.sam C.sam
#step 1 run picard
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinate
INPUT=A.sam OUTPUT=A.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinate
INPUT=C.sam OUTPUT=C.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
#step 2 mark duplicate reads
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/MarkDuplicates.jar INPUT=A.
bam OUTPUT=A.marked.bam METRICS_FILE=metrics CREATE_INDEX=true VALIDATION_
STRINGENCY=LENIENT
java -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 “-L
exonscaptured.bed” is only needed if you are doing ExomeSeq. For WGS, not
needed.
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -
R /bwa/indexes/hg19.fa -o A.bam.list -I A.marked.bam -L exonscaptured.bed
java -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 4
java -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.bed
java -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 5
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=
A.marked.realigned.bam OUTPUT=A.marked.realigned.fixed.bam SO=coordinate
VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=
C.marked.realigned.bam OUTPUT=C.marked.realigned.fixed.bam SO=coordinate
VALIDATION_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 -T
CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov
CycleCovariate -cov DinucCovariate -recalFile A.recal_data.csv -L
exonscaptured.bed
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/
hg19.fa -knownSites /DBSNP/dbsnp135.vcf -I C.marked.realigned.fixed.bam -T
CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov
CycleCovariate -cov DinucCovariate -recalFile C.recal_data.csv -L
exonscaptured.bed
#step 7
java -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.
bed
java -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 8
java -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 -L
exonscaptured.bed
#select step, select SNP
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
SelectVariants --variant A.snps.AB.vcf -o A.snpsonly.vcf -selectType SNP -L
exonscaptured.bed
#select step2, select Indel
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
SelectVariants --variant A.snps.AB.vcf -o A.indels.vcf -selectType INDEL -L
exonscaptured.bed
#step 9 more filter from HapMap and 1000 Genome project
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantRecalibrator -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 -L
exonscaptured.bed
#step10
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
ApplyRecalibration -input A.snpsonly.vcf -tranchesFile A.snp.tranches.vcf -
recalFile A.snps.recal.vcf -o A.snps.filtered.vcf -L exonscaptured.bed
#step 11
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantFiltration --variant A.indels.vcf -o A.indels.filtered.vcf --
filterExpression "QD < 2.0" --filterExpression "ReadPosRankSum < -20.0" --
filterExpression "FS > 200.0" --filterName QDFilter --filterName
ReadPosFilter --filterName FSFilter -L exonscaptured.bed
#step 12
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantFiltration --variant A.snps.filtered.vcf --clusterSize 3 --
clusterWindowSize 10 -o A.snps.finalfilteredAB.vcf -L exonscaptured.bed
我们尝试过好几种不同的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
Check this website if you have any questions:
http://seqanswers.com/wiki/How-to/exome_analysis
#step 0 align reads
bwa 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.sam
bwa 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 reads
grep chr A.sam > temp.sam
mv temp.sam A.sam
grep chr C.sam > temp2.sam
mv temp2.sam C.sam
#step 1 run picard
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinate
INPUT=A.sam OUTPUT=A.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinate
INPUT=C.sam OUTPUT=C.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
#step 2 mark duplicate reads
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/MarkDuplicates.jar INPUT=A.
bam OUTPUT=A.marked.bam METRICS_FILE=metrics CREATE_INDEX=true VALIDATION_
STRINGENCY=LENIENT
java -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 “-L
exonscaptured.bed” is only needed if you are doing ExomeSeq. For WGS, not
needed.
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -
R /bwa/indexes/hg19.fa -o A.bam.list -I A.marked.bam -L exonscaptured.bed
java -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 4
java -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.bed
java -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 5
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=
A.marked.realigned.bam OUTPUT=A.marked.realigned.fixed.bam SO=coordinate
VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=
C.marked.realigned.bam OUTPUT=C.marked.realigned.fixed.bam SO=coordinate
VALIDATION_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 -T
CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov
CycleCovariate -cov DinucCovariate -recalFile A.recal_data.csv -L
exonscaptured.bed
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/
hg19.fa -knownSites /DBSNP/dbsnp135.vcf -I C.marked.realigned.fixed.bam -T
CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov
CycleCovariate -cov DinucCovariate -recalFile C.recal_data.csv -L
exonscaptured.bed
#step 7
java -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.
bed
java -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 8
java -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 -L
exonscaptured.bed
#select step, select SNP
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
SelectVariants --variant A.snps.AB.vcf -o A.snpsonly.vcf -selectType SNP -L
exonscaptured.bed
#select step2, select Indel
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
SelectVariants --variant A.snps.AB.vcf -o A.indels.vcf -selectType INDEL -L
exonscaptured.bed
#step 9 more filter from HapMap and 1000 Genome project
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantRecalibrator -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 -L
exonscaptured.bed
#step10
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
ApplyRecalibration -input A.snpsonly.vcf -tranchesFile A.snp.tranches.vcf -
recalFile A.snps.recal.vcf -o A.snps.filtered.vcf -L exonscaptured.bed
#step 11
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantFiltration --variant A.indels.vcf -o A.indels.filtered.vcf --
filterExpression "QD < 2.0" --filterExpression "ReadPosRankSum < -20.0" --
filterExpression "FS > 200.0" --filterName QDFilter --filterName
ReadPosFilter --filterName FSFilter -L exonscaptured.bed
#step 12
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantFiltration --variant A.snps.filtered.vcf --clusterSize 3 --
clusterWindowSize 10 -o A.snps.finalfilteredAB.vcf -L exonscaptured.bed
a*7
3 楼
越狱以后可以
a*m
4 楼
现在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
，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*1
5 楼
很感激。我赶紧去对比下你的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
但问题是，貌似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
p*e
6 楼
好东西，顶一下
i*r
7 楼
谢谢，祝你一切顺利！
我这种情况怎么搞（目前还是没有弄明白）
有2个原始文件，一个是发飙了的SNPs数据，大概70个不同human poupulations，是一
个bed文件，结构是这样（右边省略了若干column）
chr1 41217 41218 snp 2 + T A dbsnp.108:
rs3863625 NN
chr1 41255 41256 snp 3 + C T dbsnp.111:
rs4543737 NN
chr1 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:0
all-hg18_2 0 chr1 138 114 9M1D86M1D22M * 0
0 （） （） AS:i:0
all-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}' | uniq
chr1
chr10
chr11
...
chr2
chr20
...
不是我想要的：
chr1
chr2
chr3
...
我这种情况怎么搞（目前还是没有弄明白）
有2个原始文件，一个是发飙了的SNPs数据，大概70个不同human poupulations，是一
个bed文件，结构是这样（右边省略了若干column）
chr1 41217 41218 snp 2 + T A dbsnp.108:
rs3863625 NN
chr1 41255 41256 snp 3 + C T dbsnp.111:
rs4543737 NN
chr1 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:0
all-hg18_2 0 chr1 138 114 9M1D86M1D22M * 0
0 （） （） AS:i:0
all-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}' | uniq
chr1
chr10
chr11
...
chr2
chr20
...
不是我想要的：
chr1
chr2
chr3
...
u*d
8 楼
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
，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*p
9 楼
攒人品，顺便回答一下 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
Check this website if you have any questions:
http://seqanswers.com/wiki/How-to/exome_analysis
#step 0 align reads
bwa 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.sam
bwa 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 reads
grep chr A.sam > temp.sam
mv temp.sam A.sam
grep chr C.sam > temp2.sam
mv temp2.sam C.sam
#step 1 run picard
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinate
INPUT=A.sam OUTPUT=A.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinate
INPUT=C.sam OUTPUT=C.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
#step 2 mark duplicate reads
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/MarkDuplicates.jar INPUT=A.
bam OUTPUT=A.marked.bam METRICS_FILE=metrics CREATE_INDEX=true VALIDATION_
STRINGENCY=LENIENT
java -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 “-L
exonscaptured.bed” is only needed if you are doing ExomeSeq. For WGS, not
needed.
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -
R /bwa/indexes/hg19.fa -o A.bam.list -I A.marked.bam -L exonscaptured.bed
java -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 4
java -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.bed
java -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 5
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=
A.marked.realigned.bam OUTPUT=A.marked.realigned.fixed.bam SO=coordinate
VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=
C.marked.realigned.bam OUTPUT=C.marked.realigned.fixed.bam SO=coordinate
VALIDATION_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 -T
CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov
CycleCovariate -cov DinucCovariate -recalFile A.recal_data.csv -L
exonscaptured.bed
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/
hg19.fa -knownSites /DBSNP/dbsnp135.vcf -I C.marked.realigned.fixed.bam -T
CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov
CycleCovariate -cov DinucCovariate -recalFile C.recal_data.csv -L
exonscaptured.bed
#step 7
java -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.
bed
java -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 8
java -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 -L
exonscaptured.bed
#select step, select SNP
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
SelectVariants --variant A.snps.AB.vcf -o A.snpsonly.vcf -selectType SNP -L
exonscaptured.bed
#select step2, select Indel
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
SelectVariants --variant A.snps.AB.vcf -o A.indels.vcf -selectType INDEL -L
exonscaptured.bed
#step 9 more filter from HapMap and 1000 Genome project
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantRecalibrator -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 -L
exonscaptured.bed
#step10
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
ApplyRecalibration -input A.snpsonly.vcf -tranchesFile A.snp.tranches.vcf -
recalFile A.snps.recal.vcf -o A.snps.filtered.vcf -L exonscaptured.bed
#step 11
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantFiltration --variant A.indels.vcf -o A.indels.filtered.vcf --
filterExpression "QD < 2.0" --filterExpression "ReadPosRankSum < -20.0" --
filterExpression "FS > 200.0" --filterName QDFilter --filterName
ReadPosFilter --filterName FSFilter -L exonscaptured.bed
#step 12
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantFiltration --variant A.snps.filtered.vcf --clusterSize 3 --
clusterWindowSize 10 -o A.snps.finalfilteredAB.vcf -L exonscaptured.bed
我们尝试过好几种不同的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
Check this website if you have any questions:
http://seqanswers.com/wiki/How-to/exome_analysis
#step 0 align reads
bwa 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.sam
bwa 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 reads
grep chr A.sam > temp.sam
mv temp.sam A.sam
grep chr C.sam > temp2.sam
mv temp2.sam C.sam
#step 1 run picard
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinate
INPUT=A.sam OUTPUT=A.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/SortSam.jar SO=coordinate
INPUT=C.sam OUTPUT=C.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
#step 2 mark duplicate reads
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/MarkDuplicates.jar INPUT=A.
bam OUTPUT=A.marked.bam METRICS_FILE=metrics CREATE_INDEX=true VALIDATION_
STRINGENCY=LENIENT
java -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 “-L
exonscaptured.bed” is only needed if you are doing ExomeSeq. For WGS, not
needed.
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -T RealignerTargetCreator -
R /bwa/indexes/hg19.fa -o A.bam.list -I A.marked.bam -L exonscaptured.bed
java -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 4
java -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.bed
java -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 5
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=
A.marked.realigned.bam OUTPUT=A.marked.realigned.fixed.bam SO=coordinate
VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
java -Xmx8g -Djava.io.tmpdir=/tmp -jar /picard/FixMateInformation.jar INPUT=
C.marked.realigned.bam OUTPUT=C.marked.realigned.fixed.bam SO=coordinate
VALIDATION_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 -T
CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov
CycleCovariate -cov DinucCovariate -recalFile A.recal_data.csv -L
exonscaptured.bed
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -l INFO -R /bwa/indexes/
hg19.fa -knownSites /DBSNP/dbsnp135.vcf -I C.marked.realigned.fixed.bam -T
CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov
CycleCovariate -cov DinucCovariate -recalFile C.recal_data.csv -L
exonscaptured.bed
#step 7
java -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.
bed
java -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 8
java -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 -L
exonscaptured.bed
#select step, select SNP
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
SelectVariants --variant A.snps.AB.vcf -o A.snpsonly.vcf -selectType SNP -L
exonscaptured.bed
#select step2, select Indel
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
SelectVariants --variant A.snps.AB.vcf -o A.indels.vcf -selectType INDEL -L
exonscaptured.bed
#step 9 more filter from HapMap and 1000 Genome project
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantRecalibrator -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 -L
exonscaptured.bed
#step10
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
ApplyRecalibration -input A.snpsonly.vcf -tranchesFile A.snp.tranches.vcf -
recalFile A.snps.recal.vcf -o A.snps.filtered.vcf -L exonscaptured.bed
#step 11
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantFiltration --variant A.indels.vcf -o A.indels.filtered.vcf --
filterExpression "QD < 2.0" --filterExpression "ReadPosRankSum < -20.0" --
filterExpression "FS > 200.0" --filterName QDFilter --filterName
ReadPosFilter --filterName FSFilter -L exonscaptured.bed
#step 12
java -Xmx8g -jar /gatk/dist/GenomeAnalysisTK.jar -R /bwa/indexes/hg19.fa -T
VariantFiltration --variant A.snps.filtered.vcf --clusterSize 3 --
clusterWindowSize 10 -o A.snps.finalfilteredAB.vcf -L exonscaptured.bed
a*m
10 楼
现在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
，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*1
11 楼
很感激。我赶紧去对比下你的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
但问题是，貌似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
p*e
12 楼
好东西，顶一下
i*r
13 楼
谢谢，祝你一切顺利！
我这种情况怎么搞（目前还是没有弄明白）
有2个原始文件，一个是发飙了的SNPs数据，大概70个不同human poupulations，是一
个bed文件，结构是这样（右边省略了若干column）
chr1 41217 41218 snp 2 + T A dbsnp.108:
rs3863625 NN
chr1 41255 41256 snp 3 + C T dbsnp.111:
rs4543737 NN
chr1 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:0
all-hg18_2 0 chr1 138 114 9M1D86M1D22M * 0
0 （） （） AS:i:0
all-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}' | uniq
chr1
chr10
chr11
...
chr2
chr20
...
不是我想要的：
chr1
chr2
chr3
...
我这种情况怎么搞（目前还是没有弄明白）
有2个原始文件，一个是发飙了的SNPs数据，大概70个不同human poupulations，是一
个bed文件，结构是这样（右边省略了若干column）
chr1 41217 41218 snp 2 + T A dbsnp.108:
rs3863625 NN
chr1 41255 41256 snp 3 + C T dbsnp.111:
rs4543737 NN
chr1 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:0
all-hg18_2 0 chr1 138 114 9M1D86M1D22M * 0
0 （） （） AS:i:0
all-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}' | uniq
chr1
chr10
chr11
...
chr2
chr20
...
不是我想要的：
chr1
chr2
chr3
...
u*d
14 楼
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
，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*r
15 楼
请问有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
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
x*d
16 楼
我用的是Tophat1.4加freebayes。 GATK老是说我的reference有问题。
u*1
17 楼
你要确定你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有问题。
reference；GATK对这个要求很严，尤其不同版本的ref，比如g1k_37和hg19；g1k会包
含一些无法assembly的contig。总之ref里的chromosome的多少，排序都很重要，同时
要注意prefix比如“chr”；20和chr20是不一样的，GATK没办法识别的
【在 x*****d 的大作中提到】
: 我用的是Tophat1.4加freebayes。 GATK老是说我的reference有问题。
x*d
18 楼
多谢指点！我用的都是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*d
19 楼
楼主少了一步吧？
bwa aln -t 4 -f -I
bwa aln -t 4 -f -I
，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
bwa aln -t 4 -f
bwa aln -t 4 -f
，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*g
20 楼
For SNP calling, there is a better one came out recently:
http://www.nature.com/ncomms/journal/v3/n12/abs/ncomms2256.html
http://www.nature.com/ncomms/journal/v3/n12/abs/ncomms2256.html
x*d
21 楼
看了,这个流程要求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
相关阅读
40岁生物千老找下家，了解行情张谢鲁事件应该可以降温了！怎么和老板处好关系？感觉老板嫌弃我嫌弃的不行了。 (转载)张鲁谢教你们怎么作人接下来就看张生家的结果能不能被重复了请问大家知道有什么比较教学和就业比较好的bioinformatics (生物信息）online Master program （全国或者湾区都行）或者湾区的part-time bioinformati生物版已死，多亏了张生天天烧纸西雅图工作机会（生物， 化学）求一本电子书若用这种羞耻的方式拿诺奖，我们中国人宁愿不要Cell signaling的抗体很难work啊翻墙求合作（丁香园）爆料！科技界最大丑闻： 饶毅和鲁白才是“谢张”事件背后的大boss！我k，张生家洗版啊这次谢张纠纷中最出格的一件事其实张谢鲁都是（中国科学）的托除了Cre/loxp跟FLP/FRT系统可以做conditional KO关于公司博后请教设计个生物过程的动画需要多少钱？张生的水军来炸版有屁用