soap2 和 soapsnp介绍

Program: SOAPaligner/soap2
Compile Date: Sun Aug 22 11:51:04 CST 2010
Author:  BGI shenzhen
Version: 2.21
Contact: soap@genomics.org.cn

Usage:  soap [options]
    -a  <str>   query a file, *.fq, *.fa
    -b  <str>   query b file
    -D  <str>   reference sequences indexing table, *.index format
    -o  <str>   output alignment file(txt)
    -M  <int>   match mode for each read or the seed part of read, which shouldn't contain more than 2 mismaches, [4]
                0: exact match only
                1: 1 mismatch match only
                2: 2 mismatch match only
                4: find the best hits
    -u  <str>   output unmapped reads file
    -t          output reads id instead reads name, [none]
    -l  <int>   align the initial n bps as a seed [256] means whole length of read
    -n  <int>   filter low-quality reads containing >n Ns before alignment, [5]
    -r  [0,1,2] how to report repeat hits, 0=none; 1=random one; 2=all, [1]
    -m  <int>   minimal insert size allowed, [400]
    -x  <int>   maximal insert size allowed, [600]
    -2  <str>   output file of unpaired alignment hits
    -v  <int>   maximum number of mismatches allowed on a read. [5] bp
    -s  <int>   minimal alignment length (for soft clip) [255] bp
    -g  <int>   one continuous gap size allowed on a read. [0] bp
    -R          for long insert size of pair end reads RF. [none](means FR pair)
    -e  <int>   will not allow gap exist inside n-bp edge of a read, default=5
    -p  <int>   number of processors to use, [1]

    -h          this help

-a  SE测序的文件或者PE测序的1个文件
-b  PE测序的另一个文件
-D ref 的索引
-o  输出文件
-M  read或者read 的部分seed 匹配的方式(不超过2个错配)
       0:只精准匹配
       1:1个mismatch
       2:   2个mismatch
       4:  找到最佳
-u  输出不匹配的reads的文件
-t  输出reads ID 而不是 name
-l  长reads的时候使用。将初始的 n-bp [256] 当做一个种子(seed) 比对,从5‘端取指定长度碱基作为种子先比对,如果能比上再比对全长的reads
-n  在比对前,过滤包含5个N 的reads 
-r  比对到多个位置的序列如何输出,0-不输出,1-随机输出一次,2-全部输出。[一般出现在repeat 区域]
-m 允许的最小插入
-x   允许的最大插入
-2  输出比对上但是不是成对reads比对上的序列的文件路径(包含文件名称)。Single-end不必输入,这个文件包含的序列都是一端被比对上的序列。
-v  一条read允许的错配数量[5]
-s  最小的比对长度(for soft clip) [255] bp
-g  一条read允许一个连续gap大小

-R     RF alignment for long insert size(>=  2k  bps)  PE  data,
                     [none] FR alignment   
        [老的解释] 插入片段大于2K设置,低于不设置

-R  for long insert size of pair end reads RF. [none](means FR pair)
      什么是 RF ? FR? 
 
-e  一条read 的n-bp 的边缘的内部不允许gap的存在
-p  线程

参数:

soap -a /workdir/SNPcall_GATK/test19610/out/524/HG2504H/HG2504H.fastq.gz -D /home/xmxjy/zsl_test/datahg19_split/hg19.fa.index -o /workdir/SNPcall_GATK/test19610/out/524/HG2504H/HG2504H.allout

输出文件格式说明:

82 TTTTCGTATGGTAAAGCCTTGGCCATTTTTGGAGCGTTTTTGGC abbbaaabbba`]aaabaaa^^a_`a`b^aZUD[aZ_^``[YO 72 a 44 +  LG_2 6082510  1 C->33G4 44M 33C10

82 TTTTCGTATGGTAAAGCCTTGGCCATTTTTGGAGCGTTTTTGGC abbbaaabbba]aaabaaa^^a_a`baZUD[aZ_``[YO 72 a 44 + LG_2 6082510 1 C->33G4 44M 33C10

输出文件:

1. 编号: read 的编号。
2. read的序列.如果read比对上参考序列的负链,会被反向互补为正链。
3. 质量值:序列的质量值,和序列顺序一致,如果read反向互补,质量值也会随着改变。
4. 比对上的次数: 最优比对的次数。没有比对上的read将被忽略。
5. a/b:pair-end比对的标记, 表示read属于来自哪个文件
6. 长度: read长度,如果是容缺失的比对,长度将是加上缺失片断的长度。
7. +/-: 比对上参考序列的正链或负链
8. 染色体名称:参考序列的染色体名称
9. 位点:第一个碱基在染色体上的位置,从1开始
10. 错配的个数
11. 错配的详细信息("C->33G4" 意思是一个错配,在参考序列的位置是第9列+33(从0开始),在参考序列上是C,read上是G,质量值是4。)
12. 比对上的数目  ("44M" 意思是44个碱基比对上了)
13. 对比的细节 ("33C10"意思是前33个比对上了,第34(参考序列上是第九列+34)个是错配,后面10个还是比对上了)

然后是soapsnp

soapsnp 官网:http://soap.genomics.org.cn/soapsnp.html

https://sourceforge.net/projects/soapsnp/

1.下载压缩包
3.解压SOAPsnp-v1.03.tar.gz 软件包
   $ tar zxvf SOAPsnp-v1.03.tar.gz 
   解压后文件夹内有两个文件SOAPsnp-v1.03.tar.gz和SOAPsnpZ.其中,SOAPsnp-v1.03.tar.gz即为原始文件,SOAPsnpZ为解压后的文件。
   (可选)为了后续操作的方便,将SOAPsnpZ更名为SOAPsnp
    $ mv SOAPsnpZ SOAPsnp
4.打开SOAPsnp,执行make all使该程序配置到系统
   $ cd SOAPsnp
   $ make all
   之后,生成的soapsnp即为安装好的程序。

1. SOAPsnp运行之前对SOAP2结果的排序

cat PE_output SE_output | sort -k8,8  ‐k9,9n > soapoutput.sort

or 
./msort -k f8 -k n9  PE_output SE_output > soapoutput.sort

msort 链接:
链接:https://pan.baidu.com/s/16ou6UH8wUMWetQNXhrDp-g
提取码:469y

参数说明:

Compulsory Parameters:
-i <FILE> Input SORTED Soap Result
-d <FILE> Reference Sequence in fasta format
-o <FILE> Output consensus file
Optional Parameters:(Default in [])
-z <Char> ASCII chracter standing for quality==0 [@]
-g <Double> Global Error Dependency Coefficient, 0.0(complete dependent)~1.0(complete independent)[0.9]
-p <Double> PCR Error Dependency Coefficient, 0.0(complete dependent)~1.0(complete independent)[0.5]
-r <Double> novel altHOM prior probability [0.0005]
-e <Double> novel HET prior probability [0.0010]
-t set transition/transversion ratio to 2:1 in prior probability
-s <FILE> Pre-formated dbSNP information
-2 specify this option will REFINE SNPs using dbSNPs information [Off]
-a <Double> Validated HET prior, if no allele frequency known [0.1]
-b <Double> Validated altHOM prior, if no allele frequency known[0.05]
-j <Double> Unvalidated HET prior, if no allele frequency known [0.02]
-k <Double> Unvalidated altHOM rate, if no allele frequency known[0.01]
-u Enable rank sum test to give HET further penalty for better accuracy. [Off]
-m Enable monoploid calling mode, this will ensure all consensus as HOM and you probably should SPECIFY higher altHOM rate. [Off]
-q Only output potential SNPs. Useful in Text output mode. [Off]
-M <FILE> Output the quality calibration matrix; the matrix can be reused with -I if you rerun the program
-I <FILE> Input previous quality calibration matrix. It cannot be used simutaneously with -M
-L <short> maximum length of read [45]
-Q <short> maximum FASTQ quality score [40]
-F <int> Output format. 0: Text; 1: GLFv2; 2: GPFv2.[0]
-E <String> Extra headers EXCEPT CHROMOSOME FIELD specified in GLFv2 output. Format is "TypeName1:DataName1:TypeName2:DataName2"[]
-T <FILE> Only call consensus on regions specified in FILE. Format: ChrName\tStart\tEnd.
-h Display this help


必须参数:
-i  输入文件,要排序
-d  ref的索引文件
-o  输出文件

可选文件:
-z   代表质量值的ASCII 码字符,0[@]
-g  <双精度型>   整体误差依赖系数,0.0(完全依赖)~1.0 (完全不依赖)[0.9]
-p  <双精度型>  PCR扩增错误依赖系数,0.0(完全依赖)~1.0 (完全不依赖)[0.5]
-r   <双精度型>  新物种纯合SNP 的先验概率[0.0005]
-e  <双精度型>  新物种杂合SNP 的先验概率[0.0010]
-t  在先验概率中,设置转换与颠换的比为2 比1
-s  dbSNPs 信息
-2 指定这个参数会使用dbSNPs信息筛选SNP[off]
-a  <双精度型>  如果不知道等位基因的频率信息,优先确认为杂合类型[0.1]
-b  <双精度型>  如果不知道等位基因的频率信息,优先确认为纯合SNP[0.05]
-j  <双精度型>  如果不知道等位基因的频率信息,不会首先确认为杂合类型[0.02]
-k  <双精度型>  如果不知道等位基因的频率信息,不会首先确认纯合类型的频率[0.01]
-u  使用秩和检验计算杂合SNP 的罚分,以提高结果的正确性[off]
-m 单倍体检测模式,这样会确保一致性组装以纯合的形式进行,并且可能需要指定高的纯合率[off]
     -n 使用二项式概率检验来确保杂合SNP 类型的准确性[off]
-q 只输出潜在的SNPs。在文本输出模式中使用。[off]
-M <文件型>  输出质量校正矩阵;如果重新运行程序,-I 参数可用来生成新的矩阵。
-I<文件型>  输入已有的质量校正矩阵。该参数不可以与-M 参数同时使用。
-L <短整型> read的最大长度[45] 。请注意,一旦有些 reads 的长度超过这个参数所设的长度,将导致程序运行的失败。
-Q <短整型> FASTQ文件的最大质量值得分[40] 。
-F <整数型>  输出格式。0:文本型;1:GLFv2 ;2:GPFv2 。[0]
-E <字符串类型>  在GLFv2 格式类型的输出时,标注被指定染色体区域的其他名字。格式如:“ 类型名字1:数据名字1:类型名字2:数据名字2”[]
-T <文件型>  只检测文件中指定范围的一致性序列。这类文件的格式为:
        染色体名字t 起始位点t 终止位点
       ……
-h 显示以上帮助 

还有


SOAPsnp 输出文件:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
chr8 35782 A R 1 A 27 1 2 G 26 1 2 5 0.500000 2.00000 1 5
chr1 1447325 G R 7 G 4 15 15 A 0 14 14 34 1.00000 1.00000 0
chr1 2786145 G R 64 A 2 66 66 G 2 22 22 132 1.00000 1.00000 0
1) Chromosome ID
2) Coordinate on chromosome, start from 1
3) Reference genotype
4) Consensus genotype
5) Quality score of consensus genotype
6) Best base
7) Average quality score of best base
8) Count of uniquely mapped best base
9) Count of all mapped best base
10) Second best bases
11) Average quality score of second best base
12) Count of uniquely mapped second best base
13) Count of all mapped second best base
14) Sequencing depth of the site
15) Rank sum test p_value
16) Average copy number of nearby region
17) Whether the site is a dbSNP.


chr8  35782  A  R  1  A  27  1  2  G  26  1  2  5   0.500000  2.00000 1  5
格式说明(从左到右)
1.   染色体名称。
2.   位点坐标。
3.   参考序列对应位置的碱基。
4.   测序个体的碱基(多种碱基会被写成简并碱基)。
5.   质量值。
6.   似然值最大的碱基(第一碱基)。
7.   似然值最大碱基的平均质量值。
8.   支持该SNP 的read 数  (只包括唯一对比上的reads)。
9.   支持该SNP 的read 数(所有 reads)。
10. 似然值次大的碱基(第二碱基)。
11. 似然值次大的碱基的质量值。
12. 支持该SNP 的read 数  (只包括唯一对比上的reads)。
13. 支持该SNP 的read 数(所有 reads)。
14. 比对上该位置的所有read 个数。
15. 秩和检验值。
16. 该位点的拷贝数估计值。
17. dbSNP 数据库中是否有该位点,1 表示存在,0 表示不存在。
18. 与此位点最近的另一个SNP 位点的距离。 (可能没有)



2 GLFv2 和GPFv2 格式
 GLFv2(基因组似然性格式v2)是一个二进制的文件,通常在群体检测 SNP
时使用。

参考:http://blog.sina.com.cn/s/blog_13de3725c0102v7pq.html

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 205,132评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,802评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,566评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,858评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,867评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,695评论 1 282
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,064评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,705评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 42,915评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,677评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,796评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,432评论 4 322
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,041评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,992评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,223评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,185评论 2 352
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,535评论 2 343

推荐阅读更多精彩内容