i***r
发帖数: 1035 | 1 我有2个文件，一个是所有某些群体的SNPs位点，另一个是最新测序的另一个群体
现在我要在新的测序里面，找是否和原已知的SNPs有交集，然后确认在新的测序里面，
那个位点是/不是SNP
请教什么工具可以做？ |
k******d
发帖数: 76 | 2 GATK 和samtools可以call SNPs。annovar好像可以比较是否有交集，最简单的是写一
个程序直接比较两个SNP文件 |
u*********1
发帖数: 2518 | 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
: 请教什么工具可以做？
|
i***r
发帖数: 1035 | 4 关键是还要根据 CIGAR code 校正每一个call， 因为我最后不仅需要知道有没有交集
，还需要知道那个SNP position 位点是不是 在我的新数据里面也是SNP
也可以自己写，但是有现成的tool会比较efficient。 而且我刚接触这个领域，容易写错
【在 k******d 的大作中提到】
: GATK 和samtools可以call SNPs。annovar好像可以比较是否有交集，最简单的是写一
: 个程序直接比较两个SNP文件
|
i***r
发帖数: 1035 | 5 我研究一下这个软件。
包子晚些发，谢谢咯
"
【在 u*********1 的大作中提到】
: 直接用annovar吧
: 把你那个”某些群体的SNP“为点作为database，但貌似要是vcf格式
: annotate_variation.pl -filter -dbtype vcf -vcfdbfile ”某些群体的SNP“.vcf "
: 另外一个群体的SNP” humandb/
: 其实如果你不做后续的比如寻找missense SNP的话，你完全可以写个很简单的script直
: 接比较好了。
|
u*********1
发帖数: 2518 | 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。 而且我刚接触这个领域，容易写错
|
i***r
发帖数: 1035 | 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。
包子先发一部分，后面继续帮助的会接着发，谢谢咯！！ |
n******7
发帖数: 12463 | 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
|
n******7
发帖数: 12463 | 9 另外感觉你是不是还想知道你的read有没有覆盖到某个snp？
这个你可以自己parse sam 文件，或者用samtools,bedtools之类的工具，我记得有算
coverage的功能 |
n******7
发帖数: 12463 | 10 另外感觉你是不是还想知道你的read有没有覆盖到某个snp？
这个你可以自己parse sam 文件，或者用samtools,bedtools之类的工具，我记得有算
coverage的功能 |
|
|
i***r
发帖数: 1035 | 11 我现在就在尝试用samtool和bedtool，没用过加上这两个tool的说明书都极其简单，而
且是好几个小的tools，不知具体该用哪个。。。
但是我相信这两个tool应该能够实现我大部分需求，甚至全部需求。 |
n******7
发帖数: 12463 | 12 简单到不至于
比如bedtools的manual
http://code.google.com/p/bedtools/downloads/detail?name=BEDTool
一开始不知道用哪个function倒是真的，试试就明白了
【在 i***r 的大作中提到】
: 我现在就在尝试用samtool和bedtool，没用过加上这两个tool的说明书都极其简单，而
: 且是好几个小的tools，不知具体该用哪个。。。
: 但是我相信这两个tool应该能够实现我大部分需求，甚至全部需求。
|
u*********1
发帖数: 2518 | 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来，然后在比
: 较。这个比较可以用一些标准工具,比如楼上提到的；或者自己写个简单的脚本。
:
: 略）
|
c*****g
发帖数: 66 | 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
|
i***r
发帖数: 1035 | 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.
:
: 略）
|
i***r
发帖数: 1035 | 16 啊？我真的觉得有点太简单（简单的意思是太过精简以至于看不懂。。。）？比如有个
mapBed的工具竟然在manual里面没有任何说明。。。
【在 n******7 的大作中提到】
: 简单到不至于
: 比如bedtools的manual
: http://code.google.com/p/bedtools/downloads/detail?name=BEDTool
: 一开始不知道用哪个function倒是真的，试试就明白了
|
w****w
发帖数: 521 | |
c*****g
发帖数: 66 | 18 多谢你的baozi！
应该不会有问题，你试试就知道了。大不了sort一下。
我猜想你应该没有单个个体的sequence data在B人群里。
如果是这样的话，那些SNP caller都不顶用的。你需要自己pileup了之后去看具体的序
列是不是有变化。pileup就是做你讲的一个个alignment地check制定位置的碱基，不用
你自己再写程序做了。
【在 i***r 的大作中提到】
: 我手上的bam 文件似乎排序有点问题（我转成sam之后，第一列是这样）：
: chr1
: chr10
: chr11
: ...
: chr2
: chr21
: ...
: chr3
: 会不会有问题？
|
w****w
发帖数: 521 | 19 如果20条reads包含某位点，其中19条是T,一条是C,quality都还可以，此位点genotype
就是TC?
【在 c*****g 的大作中提到】
: 多谢你的baozi！
: 应该不会有问题，你试试就知道了。大不了sort一下。
: 我猜想你应该没有单个个体的sequence data在B人群里。
: 如果是这样的话，那些SNP caller都不顶用的。你需要自己pileup了之后去看具体的序
: 列是不是有变化。pileup就是做你讲的一个个alignment地check制定位置的碱基，不用
: 你自己再写程序做了。
|
i***r
发帖数: 1035 | 20 pileup不行，我用你的方法，出来全是N。主要是我问题没有说清楚（发现要说清楚很
难）
就像楼上说的，pileup是对多个reads，找可能的snps
我是要找B文件里面的SNP位点，是否在A（序列）中可能存在，所以
先看B位点是否在A的区间内，如果是
再看那个位点上，A是什么碱基，如果和reference不一样，则考虑是SNP （population
A 的SNP） |
