【包子求助】call SNPs 有哪些工具？？# Biology - 生物学
i*r
1 楼
我有2个文件，一个是所有某些群体的SNPs位点，另一个是最新测序的另一个群体
现在我要在新的测序里面，找是否和原已知的SNPs有交集，然后确认在新的测序里面，
那个位点是/不是SNP
请教什么工具可以做？
现在我要在新的测序里面，找是否和原已知的SNPs有交集，然后确认在新的测序里面，
那个位点是/不是SNP
请教什么工具可以做？
k*d
2 楼
GATK 和samtools可以call SNPs。annovar好像可以比较是否有交集，最简单的是写一
个程序直接比较两个SNP文件
个程序直接比较两个SNP文件
u*1
3 楼
直接用annovar吧
把你那个”某些群体的SNP“为点作为database，但貌似要是vcf格式
annotate_variation.pl -filter -dbtype vcf -vcfdbfile ”某些群体的SNP“.vcf "
另外一个群体的SNP” humandb/
其实如果你不做后续的比如寻找missense SNP的话，你完全可以写个很简单的script直
接比较好了。
【在 i***r 的大作中提到】
: 我有2个文件，一个是所有某些群体的SNPs位点，另一个是最新测序的另一个群体
: 现在我要在新的测序里面，找是否和原已知的SNPs有交集，然后确认在新的测序里面，
: 那个位点是/不是SNP
: 请教什么工具可以做？
把你那个”某些群体的SNP“为点作为database，但貌似要是vcf格式
annotate_variation.pl -filter -dbtype vcf -vcfdbfile ”某些群体的SNP“.vcf "
另外一个群体的SNP” humandb/
其实如果你不做后续的比如寻找missense SNP的话，你完全可以写个很简单的script直
接比较好了。
【在 i***r 的大作中提到】
: 我有2个文件，一个是所有某些群体的SNPs位点，另一个是最新测序的另一个群体
: 现在我要在新的测序里面，找是否和原已知的SNPs有交集，然后确认在新的测序里面，
: 那个位点是/不是SNP
: 请教什么工具可以做？
i*r
4 楼
关键是还要根据 CIGAR code 校正每一个call， 因为我最后不仅需要知道有没有交集
，还需要知道那个SNP position 位点是不是 在我的新数据里面也是SNP
也可以自己写，但是有现成的tool会比较efficient。 而且我刚接触这个领域，容易写错
【在 k******d 的大作中提到】
: GATK 和samtools可以call SNPs。annovar好像可以比较是否有交集，最简单的是写一
: 个程序直接比较两个SNP文件
，还需要知道那个SNP position 位点是不是 在我的新数据里面也是SNP
也可以自己写，但是有现成的tool会比较efficient。 而且我刚接触这个领域，容易写错
【在 k******d 的大作中提到】
: GATK 和samtools可以call SNPs。annovar好像可以比较是否有交集，最简单的是写一
: 个程序直接比较两个SNP文件
i*r
5 楼
我研究一下这个软件。
包子晚些发，谢谢咯
"
【在 u*********1 的大作中提到】
: 直接用annovar吧
: 把你那个”某些群体的SNP“为点作为database，但貌似要是vcf格式
: annotate_variation.pl -filter -dbtype vcf -vcfdbfile ”某些群体的SNP“.vcf "
: 另外一个群体的SNP” humandb/
: 其实如果你不做后续的比如寻找missense SNP的话，你完全可以写个很简单的script直
: 接比较好了。
包子晚些发，谢谢咯
"
【在 u*********1 的大作中提到】
: 直接用annovar吧
: 把你那个”某些群体的SNP“为点作为database，但貌似要是vcf格式
: annotate_variation.pl -filter -dbtype vcf -vcfdbfile ”某些群体的SNP“.vcf "
: 另外一个群体的SNP” humandb/
: 其实如果你不做后续的比如寻找missense SNP的话，你完全可以写个很简单的script直
: 接比较好了。
u*1
6 楼
还需要知道那个SNP position 位点是不是 在我的新数据里面也是SNP
Don't quite understand. You mean "your new database"(我的新数据) is not SNP-
calling file? Then first use GATK/Samtools to call SNP/indel from "your new
database", then filter against your old database.
If you have no experience using GATK, and in a hurry to get results, I
strongly suggest using Samtools, which is basically just one bash command,
while GATK is monsterous algorithm. Also newest version of GATK is coming
out and all those old scripts may now be retired.
写错
【在 i***r 的大作中提到】
: 关键是还要根据 CIGAR code 校正每一个call， 因为我最后不仅需要知道有没有交集
: ，还需要知道那个SNP position 位点是不是 在我的新数据里面也是SNP
: 也可以自己写，但是有现成的tool会比较efficient。 而且我刚接触这个领域，容易写错
Don't quite understand. You mean "your new database"(我的新数据) is not SNP-
calling file? Then first use GATK/Samtools to call SNP/indel from "your new
database", then filter against your old database.
If you have no experience using GATK, and in a hurry to get results, I
strongly suggest using Samtools, which is basically just one bash command,
while GATK is monsterous algorithm. Also newest version of GATK is coming
out and all those old scripts may now be retired.
写错
【在 i***r 的大作中提到】
: 关键是还要根据 CIGAR code 校正每一个call， 因为我最后不仅需要知道有没有交集
: ，还需要知道那个SNP position 位点是不是 在我的新数据里面也是SNP
: 也可以自己写，但是有现成的tool会比较efficient。 而且我刚接触这个领域，容易写错
i*r
7 楼
举个例子：
数据B：已发表的我们发现的SNPs，大概这样（chromosome，位点，和SNPs，其余省略）
chr1 1240 *** C
chr1 1270 *** T
数据A：最新的测序数据，大概这样（chromosome，序列起点，末点，CIGAR，序列，其
余略）
chr1 1234 1279 * * * 20M5D20M AAAAACCCCCCTTTTTGGGGGAAAAACCCCCTTTTTGGGGG
任务
1)确定A序列包含了B里面的SNP（两个SNPs都在1234-1279的区间，是我要的）
2）进一步根据序列，找出对应位点的base
example 1 是在第6个（1240-1234=6）base上，那么是C
example 2 是在第36个（1270-1234=36）base上，根据CIGAR code，有5个deletion，
股序列应该是：
AAAAACCCCCCTTTTTGGGGG*****AAAAACCCCCTTTTTGGGGG
那么应该是也是T
我现在就是有数据A和B，需要有软件能够：自动比较区间，同时根据CIGAR code找出对
应的base 是什么type。
包子先发一部分，后面继续帮助的会接着发，谢谢咯！！
数据B：已发表的我们发现的SNPs，大概这样（chromosome，位点，和SNPs，其余省略）
chr1 1240 *** C
chr1 1270 *** T
数据A：最新的测序数据，大概这样（chromosome，序列起点，末点，CIGAR，序列，其
余略）
chr1 1234 1279 * * * 20M5D20M AAAAACCCCCCTTTTTGGGGGAAAAACCCCCTTTTTGGGGG
任务
1)确定A序列包含了B里面的SNP（两个SNPs都在1234-1279的区间，是我要的）
2）进一步根据序列，找出对应位点的base
example 1 是在第6个（1240-1234=6）base上，那么是C
example 2 是在第36个（1270-1234=36）base上，根据CIGAR code，有5个deletion，
股序列应该是：
AAAAACCCCCCTTTTTGGGGG*****AAAAACCCCCTTTTTGGGGG
那么应该是也是T
我现在就是有数据A和B，需要有软件能够：自动比较区间，同时根据CIGAR code找出对
应的base 是什么type。
包子先发一部分，后面继续帮助的会接着发，谢谢咯！！
n*7
8 楼
我没做过复杂的SNP分析，不负责任地随便说说
你数据A应该是SAM格式的alignment数据吧？不建议你直接一个read一个read的来分析
variance site，因为这个完全可能是sequencing/alignment的错误造成的。最直接可
靠的方法是用一些variance caller，比如samtools，先call出snp/indel来，然后在比
较。这个比较可以用一些标准工具,比如楼上提到的；或者自己写个简单的脚本。
略）
【在 i***r 的大作中提到】
: 举个例子：
: 数据B：已发表的我们发现的SNPs，大概这样（chromosome，位点，和SNPs，其余省略）
: chr1 1240 *** C
: chr1 1270 *** T
: 数据A：最新的测序数据，大概这样（chromosome，序列起点，末点，CIGAR，序列，其
: 余略）
: chr1 1234 1279 * * * 20M5D20M AAAAACCCCCCTTTTTGGGGGAAAAACCCCCTTTTTGGGGG
: 任务
: 1)确定A序列包含了B里面的SNP（两个SNPs都在1234-1279的区间，是我要的）
: 2）进一步根据序列，找出对应位点的base
你数据A应该是SAM格式的alignment数据吧？不建议你直接一个read一个read的来分析
variance site，因为这个完全可能是sequencing/alignment的错误造成的。最直接可
靠的方法是用一些variance caller，比如samtools，先call出snp/indel来，然后在比
较。这个比较可以用一些标准工具,比如楼上提到的；或者自己写个简单的脚本。
略）
【在 i***r 的大作中提到】
: 举个例子：
: 数据B：已发表的我们发现的SNPs，大概这样（chromosome，位点，和SNPs，其余省略）
: chr1 1240 *** C
: chr1 1270 *** T
: 数据A：最新的测序数据，大概这样（chromosome，序列起点，末点，CIGAR，序列，其
: 余略）
: chr1 1234 1279 * * * 20M5D20M AAAAACCCCCCTTTTTGGGGGAAAAACCCCCTTTTTGGGGG
: 任务
: 1)确定A序列包含了B里面的SNP（两个SNPs都在1234-1279的区间，是我要的）
: 2）进一步根据序列，找出对应位点的base
n*7
9 楼
另外感觉你是不是还想知道你的read有没有覆盖到某个snp？
这个你可以自己parse sam 文件，或者用samtools,bedtools之类的工具，我记得有算
coverage的功能
这个你可以自己parse sam 文件，或者用samtools,bedtools之类的工具，我记得有算
coverage的功能
n*7
10 楼
另外感觉你是不是还想知道你的read有没有覆盖到某个snp？
这个你可以自己parse sam 文件，或者用samtools,bedtools之类的工具，我记得有算
coverage的功能
这个你可以自己parse sam 文件，或者用samtools,bedtools之类的工具，我记得有算
coverage的功能
i*r
11 楼
我现在就在尝试用samtool和bedtool，没用过加上这两个tool的说明书都极其简单，而
且是好几个小的tools，不知具体该用哪个。。。
但是我相信这两个tool应该能够实现我大部分需求，甚至全部需求。
且是好几个小的tools，不知具体该用哪个。。。
但是我相信这两个tool应该能够实现我大部分需求，甚至全部需求。
n*7
12 楼
简单到不至于
比如bedtools的manual
http://code.google.com/p/bedtools/downloads/detail?name=BEDTool
一开始不知道用哪个function倒是真的，试试就明白了
【在 i***r 的大作中提到】
: 我现在就在尝试用samtool和bedtool，没用过加上这两个tool的说明书都极其简单，而
: 且是好几个小的tools，不知具体该用哪个。。。
: 但是我相信这两个tool应该能够实现我大部分需求，甚至全部需求。
比如bedtools的manual
http://code.google.com/p/bedtools/downloads/detail?name=BEDTool
一开始不知道用哪个function倒是真的，试试就明白了
【在 i***r 的大作中提到】
: 我现在就在尝试用samtool和bedtool，没用过加上这两个tool的说明书都极其简单，而
: 且是好几个小的tools，不知具体该用哪个。。。
: 但是我相信这两个tool应该能够实现我大部分需求，甚至全部需求。
u*1
13 楼
非常同意。
直接把你的sam用samtools来call SNP
请看：
http://samtools.sourceforge.net/mpileup.shtml
就是那个samtools/bcftools的两个command，得到一个vcf file
然后用annovar来对比这个vcf file和你的old database
貌似都用不到bedtools
【在 n******7 的大作中提到】
: 我没做过复杂的SNP分析，不负责任地随便说说
: 你数据A应该是SAM格式的alignment数据吧？不建议你直接一个read一个read的来分析
: variance site，因为这个完全可能是sequencing/alignment的错误造成的。最直接可
: 靠的方法是用一些variance caller，比如samtools，先call出snp/indel来，然后在比
: 较。这个比较可以用一些标准工具,比如楼上提到的；或者自己写个简单的脚本。
:
: 略）
直接把你的sam用samtools来call SNP
请看：
http://samtools.sourceforge.net/mpileup.shtml
就是那个samtools/bcftools的两个command，得到一个vcf file
然后用annovar来对比这个vcf file和你的old database
貌似都用不到bedtools
【在 n******7 的大作中提到】
: 我没做过复杂的SNP分析，不负责任地随便说说
: 你数据A应该是SAM格式的alignment数据吧？不建议你直接一个read一个read的来分析
: variance site，因为这个完全可能是sequencing/alignment的错误造成的。最直接可
: 靠的方法是用一些variance caller，比如samtools，先call出snp/indel来，然后在比
: 较。这个比较可以用一些标准工具,比如楼上提到的；或者自己写个简单的脚本。
:
: 略）
c*g
14 楼
这个事情很简单：
1. 根据A做一个bed file
例如
chr1 1239 1240
注意，bed是0起始的
2. samtools mpileup -l snpA.bed your.bam
and pipe the output to whatever you want.
directly reading line by line from a SAM (like you described) is a bad idea.
略）
【在 i***r 的大作中提到】
: 举个例子：
: 数据B：已发表的我们发现的SNPs，大概这样（chromosome，位点，和SNPs，其余省略）
: chr1 1240 *** C
: chr1 1270 *** T
: 数据A：最新的测序数据，大概这样（chromosome，序列起点，末点，CIGAR，序列，其
: 余略）
: chr1 1234 1279 * * * 20M5D20M AAAAACCCCCCTTTTTGGGGGAAAAACCCCCTTTTTGGGGG
: 任务
: 1)确定A序列包含了B里面的SNP（两个SNPs都在1234-1279的区间，是我要的）
: 2）进一步根据序列，找出对应位点的base
1. 根据A做一个bed file
例如
chr1 1239 1240
注意，bed是0起始的
2. samtools mpileup -l snpA.bed your.bam
and pipe the output to whatever you want.
directly reading line by line from a SAM (like you described) is a bad idea.
略）
【在 i***r 的大作中提到】
: 举个例子：
: 数据B：已发表的我们发现的SNPs，大概这样（chromosome，位点，和SNPs，其余省略）
: chr1 1240 *** C
: chr1 1270 *** T
: 数据A：最新的测序数据，大概这样（chromosome，序列起点，末点，CIGAR，序列，其
: 余略）
: chr1 1234 1279 * * * 20M5D20M AAAAACCCCCCTTTTTGGGGGAAAAACCCCCTTTTTGGGGG
: 任务
: 1)确定A序列包含了B里面的SNP（两个SNPs都在1234-1279的区间，是我要的）
: 2）进一步根据序列，找出对应位点的base
i*r
15 楼
我手上的bam 文件似乎排序有点问题（我转成sam之后，第一列是这样）：
chr1
chr10
chr11
...
chr2
chr21
...
chr3
会不会有问题？
idea.
【在 c*****g 的大作中提到】
: 这个事情很简单：
: 1. 根据A做一个bed file
: 例如
: chr1 1239 1240
: 注意，bed是0起始的
: 2. samtools mpileup -l snpA.bed your.bam
: and pipe the output to whatever you want.
: directly reading line by line from a SAM (like you described) is a bad idea.
:
: 略）
chr1
chr10
chr11
...
chr2
chr21
...
chr3
会不会有问题？
idea.
【在 c*****g 的大作中提到】
: 这个事情很简单：
: 1. 根据A做一个bed file
: 例如
: chr1 1239 1240
: 注意，bed是0起始的
: 2. samtools mpileup -l snpA.bed your.bam
: and pipe the output to whatever you want.
: directly reading line by line from a SAM (like you described) is a bad idea.
:
: 略）
i*r
16 楼
啊？我真的觉得有点太简单（简单的意思是太过精简以至于看不懂。。。）？比如有个
mapBed的工具竟然在manual里面没有任何说明。。。
【在 n******7 的大作中提到】
: 简单到不至于
: 比如bedtools的manual
: http://code.google.com/p/bedtools/downloads/detail?name=BEDTool
: 一开始不知道用哪个function倒是真的，试试就明白了
mapBed的工具竟然在manual里面没有任何说明。。。
【在 n******7 的大作中提到】
: 简单到不至于
: 比如bedtools的manual
: http://code.google.com/p/bedtools/downloads/detail?name=BEDTool
: 一开始不知道用哪个function倒是真的，试试就明白了
w*w
17 楼
Partek最容易，不过license一年要$6k多。
http://www.partek.com/Tutorials/microarray/NextGen/DetectingSNP
http://www.partek.com/Tutorials/microarray/NextGen/DetectingSNP
c*g
18 楼
多谢你的baozi！
应该不会有问题，你试试就知道了。大不了sort一下。
我猜想你应该没有单个个体的sequence data在B人群里。
如果是这样的话，那些SNP caller都不顶用的。你需要自己pileup了之后去看具体的序
列是不是有变化。pileup就是做你讲的一个个alignment地check制定位置的碱基，不用
你自己再写程序做了。
【在 i***r 的大作中提到】
: 我手上的bam 文件似乎排序有点问题（我转成sam之后，第一列是这样）：
: chr1
: chr10
: chr11
: ...
: chr2
: chr21
: ...
: chr3
: 会不会有问题？
应该不会有问题，你试试就知道了。大不了sort一下。
我猜想你应该没有单个个体的sequence data在B人群里。
如果是这样的话，那些SNP caller都不顶用的。你需要自己pileup了之后去看具体的序
列是不是有变化。pileup就是做你讲的一个个alignment地check制定位置的碱基，不用
你自己再写程序做了。
【在 i***r 的大作中提到】
: 我手上的bam 文件似乎排序有点问题（我转成sam之后，第一列是这样）：
: chr1
: chr10
: chr11
: ...
: chr2
: chr21
: ...
: chr3
: 会不会有问题？
w*w
19 楼
如果20条reads包含某位点，其中19条是T,一条是C,quality都还可以，此位点genotype
就是TC?
【在 c*****g 的大作中提到】
: 多谢你的baozi！
: 应该不会有问题，你试试就知道了。大不了sort一下。
: 我猜想你应该没有单个个体的sequence data在B人群里。
: 如果是这样的话，那些SNP caller都不顶用的。你需要自己pileup了之后去看具体的序
: 列是不是有变化。pileup就是做你讲的一个个alignment地check制定位置的碱基，不用
: 你自己再写程序做了。
就是TC?
【在 c*****g 的大作中提到】
: 多谢你的baozi！
: 应该不会有问题，你试试就知道了。大不了sort一下。
: 我猜想你应该没有单个个体的sequence data在B人群里。
: 如果是这样的话，那些SNP caller都不顶用的。你需要自己pileup了之后去看具体的序
: 列是不是有变化。pileup就是做你讲的一个个alignment地check制定位置的碱基，不用
: 你自己再写程序做了。
i*r
20 楼
pileup不行，我用你的方法，出来全是N。主要是我问题没有说清楚（发现要说清楚很
难）
就像楼上说的，pileup是对多个reads，找可能的snps
我是要找B文件里面的SNP位点，是否在A（序列）中可能存在，所以
先看B位点是否在A的区间内，如果是
再看那个位点上，A是什么碱基，如果和reference不一样，则考虑是SNP （population
A 的SNP）
难）
就像楼上说的，pileup是对多个reads，找可能的snps
我是要找B文件里面的SNP位点，是否在A（序列）中可能存在，所以
先看B位点是否在A的区间内，如果是
再看那个位点上，A是什么碱基，如果和reference不一样，则考虑是SNP （population
A 的SNP）
相关阅读
Nature再发劝退文吐槽一个bioinformatics老板很庆幸，你没有转行 (转载)FDA批准抗癌新药联系博士后的问题，前辈请进求文献一篇，没有包子只有真诚祝福！After get extraction, why DNA quality is very bad? Can I still use it?中科院苏州医工所显微光学方向招聘 (转载)求Protein & Peptide Letters文献一篇promoter的甲基化或者乙酰化可不可以用Luc reporter来检测SAS Programmer Opportunity求推荐做lentivirus的公司deregulated pathway版上有木有牛人PhD阶段 就开始审稿的？生物医学从业人员等级分类中秋快乐！paper help大家觉得线虫领域的rising star有哪些啊上位遗传突变分析高手请进：这种情况肿么破求paper 多谢。