融合基因

tophat-fusion软件分析说明

TopHat-Fusion可以从单端或双端的RNA-seq数据中识别融合基因。它将首次未比对到基因组上的reads进行打断,通过重新比对寻找候选的融合基因,后续采用多重规则进行过滤和筛选。该算法不依赖于已知的基因注释信息,这使其具有更高的敏感性,除了鉴定已知基因的融合,它还可以发现来自新基因或已知基因新剪切体的融合产物

TopHat-Fusion下载

下载 Bowtie indexes (Bowtie1或者 Bowtie2), refGene.txt, and ensGene.txt .

从 blast database中下载 human_genomic, other_genomic, and nt* .

下载数据库(blastdb,ftp://ftp.ncbi.nlm.nih.gov/blast/db/)

使用TopHat-Fusion
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#####一个例子
tophat -o tophat_MCF7 -p 8 --fusion-search --keep-fasta-order --bowtie1 --no-coverage-search -r 0 \
--mate-std-dev 80 --max-intron-length 100000 --fusion-min-dist 100000 --fusion-anchor-length 13 \
--fusion-ignore-chromosomes chrM /path/to/h_sapiens/bowtie_index SRR064286_1.fastq SRR064286_2.fastq

#####参数说明

-o:用tophat_(sample_name)作为输出目录
-p:线程数
--fusion-search 打开融合基因检测命令
--bowtie1 调用bowtie1
--coverage-search 关闭--coverage-search因为它会占用大量内存且令程序运行缓慢
-r inner mate distance
--mate-std-dev inner mate distance 的标准偏差,TopHat-Fusion利用这部分区域检测融合基因
--fusion-min-dist 由于染色体重排产生的染色体内部融合最小距离
-fusion-anchor-length reads比对锚定长度
#####在tophat运行之后, 运行 tophat-fusion-post 来筛选出融合候选基因

tophat-fusion-post -p 8 --num-fusion-reads 1 --num-fusion-pairs 2 --num-fusion-both 5 /path/to/h_sapiens/bowtie_index

#####tophat-fusion-post可选参数

-v/--version 打印出帮助信息和退出.
-p/--num-threads 线程数,默认值为1.
--num-fusion-reads 支持融合基因事件最少的reads数,默认值为3.
--num-fusion-pairs 支持融合基因事件最少的双融合数,默认值为2.
--num-fusion-both 支持融合基因事件最少的双融合数和reads数的总数,默认值为0.
--max-num-fusions 最大的融合基因数目,默认值为500.
--fusion-read-mismatches mismatch值,默认值为2.
--fusion-multireads 忽略的具有多比对位点的Reads 默认值为2.
--non-human 注释的为非人的项目使用此参数.
--fusion-pair-dist 特定的融合距离,默认值为250.

SOAPfuse分析软件说明

SOAPfuse算法首先通过比对到基因组和转录本中双末端(pair end)关系的序列寻找候选的基因融合,然后采用局部穷举算法和一系列精细的过滤策略,在尽量保留真实融合的情况下过滤掉其中假阳性的基因融合。同时,该算法还具有融合断点预测和可视化功能,这对临床分子分型和肿瘤新药的开发具有重要意义。

软件安装

1.下载软件包http://soap.genomics.org.cn/soapfuse.html

1
2
$ tar -xzf /PATH_WHERE_YOU_PUT_THE_TARBALL/SOAPfuse-vX.X.tar.gz
$ cd SOAPfuse-vX.X/

Prepare

1.准备样品列表

基于下面的格式列表文件准备样品(四列)。

1 2 3 4
[sample_ID] [sequence_library_ID] [run_ID] [read_length]

2.准备RNA-Seq数据

根据样品列表存储fastq/fasta格式文件

构建目录存储RNA-Seq数据的五个要求:

(1)在主目录下存储所有RNA-Seq数据文件的子目录。我们称之为“WHOLE_SEQ-DATA_DIR”。我们公司这步的文件夹命名为“seq_data”

(2)使用sample_ID 命名 WHOLE_SEQ-DATA_DIR下的子目录的名称。我们称之为“SAMPLE_DIR”,我们公司这步的文件夹命名与样品名一致。

(3)使用sequence_library_ID 命名SAMPLE_DIR的子目录的名称。我们称之为“LIB_DIR”,与公司一致。

(4)RNA测序数据(FASTQ/ FASTA)文件存储在LIB_DIR中以其Run_ID作为文件名前缀,与公司一致。
SOAPfuse处理paired-end reads,所以前缀也应该含有双端reads序列号,如“Run_ID_1”、“Run_ID_2”。

(5)对 RNA-Seq 数据文件后缀没有要求,一般我们使用’fq.gz’ (fastq) 或者 ‘fa.gz’ (fasta)。

Run SOAPfuse

SOAPfuse运行,需要准备配置文件,SOAPfuse将根据配置文件来运行。

1.检查配置文件:

1
2
$ cd /PATH_WHERE_YOU_PUT_THE_PACKAGE/SOAPfuse-vX.X/config/,我们公司路径为/PUBLIC/software/HUMAN/SOAPfuse-v1.26/config
$ less -S config.txt

注意:

(1)所有以“#”开头的行都被视为注释。

(2)值和参数名称由“=”隔开,’=’后为值,可以被修改。

(3)根据需要某些值可以设置为“是”或“否”,有些可以保留为默认值。

(4)检查每个参数的前缀

有五种前缀,分别是“DB”,“PG”,“PS”,“PD”和“PA”。

a.“DB”指数据库的信息。

b.“PG”是指程序的信息。

c.“PS”是指流程步骤的信息。

d.“PD”是指流程目录的信息。

e.“PA”是指参数的信息。

#“DB”,“PG”,“PS”和“PD”与数据库相关,所以一旦这些参数准确地设置SOAPfuse能成功地运行。 “PA”与每个步骤的参数相关,并已设置为默认值,这样你就可以在你的第一次尝试时忽略它们,选择默认值。但是,“PA_all_fq_postfix”,它定义了RNA测序数据文件后缀,在运行前应该根据自己的RNA-Seq 文件设置参数值。

2.Run SOAPfuse

1
2
3
4
5
6
7
8
9
10
perl SOAPfuse-RUN.pl -c <config_file> -fd <WHOLE_SEQ-DATA_DIR> -l <sample_list> -o <out_directory> [Options]

#####参数说明
-c 运行流程运用的配制文件. <required>
-fd 存储双端测序Reads文件的目录 . <required> Sequenced Reads 可以为fastq 或者 fasta格式文件
-l 要处理的样本列表信息. <required>
-o 结果存储目录
-fs 想要开始的步骤. [1]
-es 想要结束的步骤. [9] 第9步是 SOAPfuse pipeline的最后一步
-h 显示帮助信息.

添加lncRNA癌基因数据库注释

再lncRNA的差异分析中添加癌基因注释,本次使用lncRNA癌基因数据库为http://www.cuilab.cn/lncrnadisease

注释脚本如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#!/usr/bin/env python
#coding=utf-8

import argparse
parser=argparse.ArgumentParser(prog='annotation',usage='%(prog)s [opthions] [value]',description='This program is used to annotate the differential expression genes!')

parser.add_argument('-D','--diffgene',help='the diffgene file')
parser.add_argument('-B','--database',help='the annotation file')
parser.add_argument('-O','--outfile',help='the output file')
argv=vars(parser.parse_args())

if argv['diffgene'] == None:
raise Exception('You should provide the diffgene file!')
if argv['diffgene'] != None:
diffgene=argv['diffgene'].strip()

if argv['database'] == None:
raise Exception('You should provide the annotation file!')
if argv['database'] != None:
database=argv['database'].strip()

if argv['outfile'] == None:
raise Exception('You should provide the output file!')
if argv['outfile'] != None:
outfile=argv['outfile'].strip()

genelist={}
Diffgene=open(diffgene)
header=Diffgene.readline()
for eachline in Diffgene:
eachlines=eachline.strip().split('\t')
genelist[eachlines[1]] = eachlines[0]
Diffgene.close()


Annotation=open(database)
header=Annotation.readline()
Outfile=open(outfile,'w')
Outfile.write(header)

for eachline in Annotation:
eachlines=eachline.strip().split('\t')
if eachlines[0] in genelist.keys():
Outfile.write(genelist[eachlines[0]] + "\t" + eachline)

Annotation.close()
Outfile.close()