由买买提看人间百态

boards

本页内容为未名空间相应帖子的节选和存档，一周内的贴子最多显示50字，超过一周显示500字 访问原贴
Biology版 - 贡献一个SNP/Indel calling pipeline
相关主题
● 版上有谁用过或知道Knome这个公司吗?● 做线虫的有人试过WGS+SNP的方法一步测序出mutant吗？
● bioinformatics吐下槽● machine learning来对GWAS结果建模
● 有谁谈谈从零开始学NGS数据分析都需要具备什么知识？● 问个whole exome capture之后出来的data要怎么分析
● 该转到computational bio领域吗● 请教Bioinformatics职业规划~~~
● 新手请教CNV caller● bioinformatics postdoc poition($35,000 - $40,000)
● NGS_Illumina类● 急问：需要多少内存
● 【包子求助】call SNPs 有哪些工具？？● Bioinformatics招人 提供refer
● 请教个DNA相关的实验问题● 一半实验一半生物信息求选择建议
相关话题的讨论汇总
话题: gatk话题: bwa话题: xmx8g话题: jar话题: indexes
进入Biology版参与讨论
1 (共1页)
j*p
发帖数: 411
1
攒人品，顺便回答一下 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
a**m
发帖数: 184
2
现在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*********1
发帖数: 2518
3
很感激。我赶紧去对比下你的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
p**********e
发帖数: 20
4
好东西，顶一下
i***r
发帖数: 1035
5
谢谢，祝你一切顺利！
我这种情况怎么搞（目前还是没有弄明白）
有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
发帖数: 573
6
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*p
发帖数: 411
7
攒人品，顺便回答一下 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
a**m
发帖数: 184
8
现在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*********1
发帖数: 2518
9
很感激。我赶紧去对比下你的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
p**********e
发帖数: 20
10
好东西，顶一下
相关主题
● NGS_Illumina类● 做线虫的有人试过WGS+SNP的方法一步测序出mutant吗？
● 【包子求助】call SNPs 有哪些工具？？● machine learning来对GWAS结果建模
● 请教个DNA相关的实验问题● 问个whole exome capture之后出来的data要怎么分析
进入Biology版参与讨论
i***r
发帖数: 1035
11
谢谢，祝你一切顺利！
我这种情况怎么搞（目前还是没有弄明白）
有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
发帖数: 573
12
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******r
发帖数: 91
13
请问有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
x*****d
发帖数: 704
14
我用的是Tophat1.4加freebayes。 GATK老是说我的reference有问题。
u*********1
发帖数: 2518
15
你要确定你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*****d
发帖数: 704
16

多谢指点！我用的都是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
发帖数: 704
17
楼主少了一步吧？
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*****g
发帖数: 78
18
For SNP calling, there is a better one came out recently:
http://www.nature.com/ncomms/journal/v3/n12/abs/ncomms2256.html
x*****d
发帖数: 704
19

看了,这个流程要求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
1 (共1页)
进入Biology版参与讨论
相关主题
● 一半实验一半生物信息求选择建议● 新手请教CNV caller
● 制药公司招生物信息Senior Information Scientist● NGS_Illumina类
● 高年级PhD毕业求建议● 【包子求助】call SNPs 有哪些工具？？
● 下一代技术测序分析结果需要会什么软件技术？● 请教个DNA相关的实验问题
● 版上有谁用过或知道Knome这个公司吗?● 做线虫的有人试过WGS+SNP的方法一步测序出mutant吗？
● bioinformatics吐下槽● machine learning来对GWAS结果建模
● 有谁谈谈从零开始学NGS数据分析都需要具备什么知识？● 问个whole exome capture之后出来的data要怎么分析
● 该转到computational bio领域吗● 请教Bioinformatics职业规划~~~
相关话题的讨论汇总
话题: gatk话题: bwa话题: xmx8g话题: jar话题: indexes

未名新帖统计// 7月16日

#版面帖数(主题数)
-全站4871 (796)
1Military3777 (569)
2Stock341 (51)
3Joke117 (17)
4History116 (3)
5Automobile100 (9)
6USANews55 (9)
7Midlife45 (1)
8Headline41 (41)
9Dreamer33 (13)
10FleaMarket32 (20)
11Living30 (7)

* 这里只显示发帖超过25的版面，努力灌水吧:-)

历史上的今天

  1. faintcat妹妹看进来～～ 发表于12年前.
  2. NSC, PD 1/7/2007, EB2, ... 发表于11年前.
  3. [FBA求购]MJVE2 758 MJVM2 ... 发表于6年前.
  4. 老生常谈,归与不归 发表于10年前.
  5. 【申请】Seattle西雅图 版版主——申请人... 发表于9年前.
  6. 宝宝出生，头骨骨折，求祝福 发表于9年前.
  7. 求推荐舒缓优美的古典音乐 发表于11年前.
  8. 百分之一的北京人上北大 中国网友愤怒(转载） 发表于10年前.
  9. 新人带狗狗Bailey来报道 发表于12年前.
  10. 全世界最有价值的运动队 发表于10年前.
  11. 请问大切诺基的质量如何 发表于6年前.
  12. TNND，军版全是BKC 发表于15年前.
  13. Inception 发表于12年前.
  14. 微软的有些家属可真恶心，为了卖保险脸都不要了 发表于10年前.
  15. 每周坐高铁的苦逼来说说感受吧!! 发表于9年前.
除非另有声明，本站内容采用Creative Commons BY-NC-SA 3.0协议进行许可，转载请注明来自未名观察 - 隐私政策2011-07-24 10:06:12由admin编辑