一、 运行 meerkat
前面已经依序 meerkat 的环境和 meerkat,运行了预处理一步,在相对应的 bam 文件目录下生成了大批文件,因此,当要用 meerkat 处理某个 bam 文件时,应先将该 bam 文件移动到专有的一个文件夹,manual 中也建议这样用。
预处理生成的文件包括:
黑名单文件. gz
isinfo 文件:包括插入大小信息
pdf 文件:插入大小的分布图,unmapped reads 长度的分布图,softclip reads 长度分布图
pre.log 文件:日志文件,包括输入的参数,输入输出信息,reads 数,unallign reads 数等
softclips.fq.gz: softlip reads 文件,完整的 read
softclips.rdist:softclip 的 reads 长度分布信息记录
unmapped.fq.gz:unmapped 的 reads 文件,完整的 read
unmapped.rdist: unmapped 的 reads 长度的分布信息
sr1.fq.gz : 从 softclip read 或者 unmaped read 切下来的指定 bp 的 reads
sr2.fq.gz : 与 sr1.fq 配对的 reads
一个文件夹包括 1_1fq.gz,1_2fq.gz , 里面有不同长度的 reads。每个 read group 人工的 reads 对
1.1 meerkat.pl
perl ./scripts/meerkat.pl [options]
-b FIle 已经 sort 过的 bam 文件
-k int [0/1] 是否使用预处理产生的黑名单文件,默认 1
-d FLT call discordant mate pairs 的标准差阈值,默认 3. 这个选项控制怎样寻找 discordant read 对,等同于定义最大的 concordant fragment(配对的 reads 定位到的片段),计算方法是:中位插入大小 + d x sd, 如果插入大小分布图狭窄并对称,就使用默认参数,例如下面一二图。如果分布图偏宽,或者带着长尾,三四图,参数应设为 5,对于高深度 (>30x),尽管分布图狭窄,但是使用 5 会轻微好一点。如果峰窄,但是仍然带着一个尾部 (图五),好像大部分 TCGA 数据那样,在 meekat.pl 这一步使用 5 比 3 更好
-c FLT 集簇 discordant mate pair 标准差阈值,默认和 - d 相同。控制怎样集簇,构建断点的置信区间。等同于定义覆盖断点的最大片段(中位插入大小 + c x sd)。如果 - d 选项是 5 或者小于 5,使用用 - c 相同的参数,如果 - d 很大,比如 10,那么使用更小的 - c,比如 5。如果数据有很高的测序深度,或者有广泛的拷贝的变化,那么使用 5 而不是 3 来避免不能构建断点的置信区间。
-p FLT call 一个事件支持的配对数阈值,默认 2
-o INT 需要支持全长 reads 对的数目,默认是 0,设定这个选项会降低小复杂事件的敏感度。如果使用 - a 1 -l 1 推荐使用 - o 1 来避免重复序列导致的过多人工产物
-q INT call 一个事件支持的 split reads 数,默认是 1
-z INT 事件数的最大值,默认 1,000,000,000
-s INT 从 unmapped reads 开始和末端切下来的 bp 数,必须和与处理的 - s 参数设置一样
-m INT [0/1],如果设定为 1,使用 meerkat 去去除 duplicates,如果设为 0,使用 Picard 的 flag 'd' marked ,或者其他工具去除 duplicates.bwa mem 的数据必须用 picard mark duplicate ,所以要设为 0.
-u INT [0/1], 使用 bam 文件的全部比对,如果 BAM 文件不是由 BWA 产生的,或者由 bwa 产生但是没有 XT 标签,那么开这个选项,开了这个选项会强制关闭 - a 选项。默认是 0。开了这个选项之后应使用 - Q 参数过滤低 mapping reads, 推荐 - Q 设置 10,对于 bwa 比对过的 bam,也可以继续过滤。
-Q INT 被使用的 reads 的最小 mapping quality,默认是 0
-g INT 在主要的 bam 文件使用的备择 mapping 数,默认使用全部。备择 mapping 数之前通过 bwa -N 参数输出到 XA 标签中。bwa mem 默认 输出备择 mapping。
-f INT 在 clipped 比对中考虑的备择 mapping 数,默认使用全部。
-l INT [0/1], 是否考虑 clipped 比对,默认 1. 和预处理一样,对于 bwa mem 出来的基因组,不需要重新 mapping.
-t INT 在 bwa 比对中用到的核数,默认 1
-R STR 包括黑名单 read group 的文件,每一行一个 read group ID
-F STR fasta 文件夹路径
-S STR samtools 路径,如果 samtools 不在环境变量中的话
-W STR bwa 路径,如果 bwa 不在环境变量中的话
-B STR blastall 和 formatdb 的路径,如果不在环境变量的话
-P STR 指定运行顺序,dc|cl|mpd|alg|srd|rf|all,默认 all,每一步运行都需要上一步的结果
dc: extract discordant read pairs,提取 discordant read 对
cl: construct clusters of discordant read pairs,将 discordant read 对建簇
mpd: call events based on read pairs,基于 read 对 call 事件
alg: align split reads to candidate break point regions,比对 split read 到候选的断点区域。如果你有多核心的话,可以将 alg 步骤切为两步, alg1 和 alg2,alg1 运行多线程,alg2 运行单线程。
srd: confirm events based on split reads and filter results,基于 split reads 和过滤的结果对事件进行证明
rf: refine break points by local alignments,通过区域比对定义断点
all: run all above steps , 运行全部步骤
-h 帮助
manual 中的举例:
50bp reads,<10x 的 TCGA 基因组 使用 - s 18 -d 5 -a 0 -l 0 -q 1,猜测:reads 长度较小,所以取 1/3 read 长度,-s 取 18, TCGA 基因组, 插入分布狭窄带尾,所以 - d 设为 5, 测序深度较低,reads 长度较短,所以尽量保留数据,-a 设为 0, -l 设为 0,-q 设的较低,1。
101bp reads, 60-80x TCGA 基因组 使用 - s 20 -d 5 -p 5 -o 3,75-101bp 使用 - s 20,TCGA 基因组使用 - d 5,测序深度深,-o 设置更高 3。
如果肿瘤基因组 60x,相对应的正常基因组 30x,那么使用 60x 的参数,常常用配对的正常组织中的 black.list.gz 重命名并放到肿瘤 bam 文件处理的文件夹,替换 cancer 的 blacklist.gz 文件。
1.2 mechanism.pl
perl ./scripts/mechanism.pl [options]
-b FILE sort 过的 bam 文件
-o INT [0/1] 在 TE 包括 rmsk 类型 \"Other\",默认 1
-t INT TE 的最大值,默认 100000
-z INT 被处理的 SV 的数量限制,默认 1000000000
-R STR repeat 注释路径,能够从 UCSC 下载
-h help
二、参照
manual 中给出的数据,HapMap 个体 NA18507(42x 深度,100bp 读长,500bp 插入大小),用 10 核 bwa 比对花费 1.5 天和 10GB 的内存。30x 深度的肿瘤数据,要多于两天并且超过 30G 的内存,如果测序质量不太好,比如很多的嵌合比对和许多非单一 mapped 的 reads,就会花费更多的时间和内存。、
三、结果
运行完预处理和 meerkat.pl 后,能够得到两个文件. intra.refined.typ.sorted 和 inter.refined.typ.sorted,运行完 mechanism.pl 后,会得到. variants 文件,否则,就该回去看一下设置是否出现问题。
运行的肿瘤基因组文件的时候,一定要把正常组织的 blacklist 文件替换肿瘤基因组的 blacklist 文件,blacklist 文件可以在预处理中生成。如果在预处理中出现错误信息 "differing read lengths",没有关系,仅仅是告诉你在一些 read group 中长度不一致。
四、包含的其他程序
4.1 转换 variant 文件到 vcf 格式
- perl ./scripts/meerkat2vcf.pl -i variantfile -o vcffile -H headerfile -F /db/hg18/hg18_fasta/
-F 选项可以产生一个头文件
4.2 融合基因注释
- perl ./scripts/fusions.pl -i variantfile -G /db/hg18/refGene_hg18_sorted.txt
4.3 为变异位点设计引物
使用 primer.pl
-i FILE 输入文件,fusion.pl 产生的融合文件
-o FILE 输出文件
-p STR 引物固定序列
-c INT 列数补偿,默认 0,如果第一列存在如何事件的样品名称,那么设为 1
-f INT 侧翼区域,默认 500,设计引物所在的区域
-s STR 被 "," 分隔的引物大小,默认 20,23,25,27
-m STR 引物最小,最优,最大 Tm 值,"," 分隔,默认 50,60,65
-n INT 对每一个引物片段,设计的音物个数,默认 5
-r INT 掩盖 repeat,默认 0
-q INT 输出设计引物的侧翼序列
-F STR fasta 文件文件夹路径
-P STR primer_core 程序文件夹路径
-L STR bla 文件夹 t 路径
-V STR blat 服务器,例如 fServer start 10.11.240.76 17777/reference/hg18/hg18.2bit -stepSize=5,服务器路径应该为 10.11.240.76
-T STR blat 端口,上面例子中的 17777
-h help
全部音物都是由 Primer3 生成,对于每一个事件,挑出. 1 和. 2,不同的取向认为是不同的事件,所以取引物的时候直接拷贝出来,不需要额外的反向互补,如果序列是小写字母,那么说明引物在重复序列。理论上,你应该用 1 blat hit 挑出小写的引物。如果 blat hit 是 0,意味着由太多 hit 了,所以不要挑选这种引物。有时候,即使引物是在重复序列上(小写字母),但是在基因组上仍然是单一比对的,(1 blat hit),因为是重复元件的变异,挑选这种引物是可以的。如果由一些引物 PCR 没有结果,你可以挑选 2 个正向引物,两个反向引物来同时进行 4 个 PCR 反应。引物设计的普遍规则仍然要使用,比如,你应该挑选 TM 值相差不大并且 GC 含量不太极端的。
4.4 计算等位基因频率
使用 discon.pl, 这个脚本会告诉你全部断点的 discordant 和 concordant 的 read 对.
-i FILE 输入 variants 文件,必须
-o FILE 输入文件,必须
-D INT 从 bam 文件中计算 discordant pair 的数目,默认 0,基因分型的时候打开这个选项
-B FILE bam 文件,必须
-C FILE 由 Meerkat 产生的 cluster 文件,必须
-I FILE Meerkat 运行时产生的 isinfo 文件
-K FILE 包含要忽略的 read group 的文件名,一个 read ID 一行,和 meerkat.pl 的 R 参数一样
-S FILE samtools 文件夹路径,如果 samtools 不再环境变量中
-d FLT call discordant read 对的标准差阈值,默认 3
-Q INT 使用的 reads 最小的 mapping quality,默认 0
-h help
- perl ./scripts/discon.pl -d 5 -Q 10 -i $somaticg -o $somaticg_rp -B $cancer_bam -C
- $cancer_cluster -I $cancer_isinfo -K $cancer_blacklistrg -S /home/ly55/opt/samtools/
结果文件里面每一个事件有一个 RP tag ,A_B_C_D 格式
A:全长的 discordant read 对数
B:从部分比对的 reads 中 discordant 的 reads 数
C:第一断点的 concordant reads 数
D:第二断点的 concordant reads 数
A+B 应该等与 Meerkat 得到的总的 reads 数
以上内容仅作个人笔记使用,仅供参考
参考资料
meerkat manual :http://gensoft.pasteur.fr/docs/Meerkat/0.189/Manual_0.189.pdf
来源: http://www.bubuko.com/infodetail-1974354.html