转录本融合位点上下游序列获取

转录本融合位点上下游序列获取
展开全文
首先需要理解：基因，转录本
(transcripts,isoform，mRNA序列)、EXON
区域，cDNA序列、UTR区域，ORF序列、
CDS序列这些概念，一个基因可以转录为多
个转录本，真核生物里面每个转录本通常是由
一个或者多个EXON组成，能翻译为蛋白的
EXON区域是CDS区域，不能翻译的那些
EXON的开头和结尾是UTR区域，翻译区域
合起来是ORF序列，而转录本逆转录就是
cDNA序列。

STAR-Fusion 结果解读
通常建议大家对RNA-seq数据使用 STAR-Fusion 来检测转录本融合现象，得到的结果如下：
文件'star-fusion.fusion_predictions.abridged.tsv', with the following format:
#FusionName JunctionReadCount SpanningFragCou nt SpliceType LeftGene LeftBreakpoint Ri ghtGene RightBreakpoint LargeAnchorSupport FFPM LeftBreakDinuc LeftBreakEntropy RightBreakDinuc RightBreakEntropy annots
THRA--
AC090627.1 27 93 ONLY_REF_SPLICE T HRA^ENSG00000126351.8chr17:38243106:+ AC090627.1 ^ENSG00000235300.3chr17:46371709:+ YES_LDAS 2 3875.8456GT 1.8892AG 1.9656 ["CCLE","FA_CancerSupp","INTRACHROMOSOMAL[chr17:8.12M
b]"]
可以看到列比较多，值得我们关心的是转录本融合的左右两个转录本的ID及其融合位点在基因组的坐标，如下：
THRA^ENSG00000126351.8chr17:38243106:+
AC090627.1^ENSG00000235300.3chr17:46371709:+
上面的坐标可以无限多，文件命名为pos.txt吧，后面写脚本需要用到，值得注意的是作者软件举例应该是hg19的基因组坐标，目前几乎都是hg38了，所以后面我脚本也是基于hg38的。

这个时候我们可能需要设计引物来对该融合转录本进行验证，所以会需要这个融合点左右两个基因的指定转录本的cDNA序列。

我们很容易拿到各个转录本的基因组坐标，但是融合点的基因组坐标不能简单对应到转录本cDNA序列里面坐标。

我们的突破点，就是找到融合点的基因组坐标到底对应到转录本cDNA序列的哪个位置。

物种的基因及其注释文件的下载
这里首选：/info/data/ftp/index.html
wget ftp:///pub/release-
95/fasta/homo_sapiens/cds/Homo_sapiens.GRCh38.cds.all.fa.gz wget ftp:///pub/release-
95/gtf/homo_sapiens/Homo_sapiens.GRCh38.95.chr.gtf.gz
首先解析基因组注释文件，就是gtf的：
1 havana exon 11869 12227 . + . gene_i
d "ENSG00000223972"; gene_version "5"; transcript_id "ENST00 000456328"; transcript_version "2"; exon_number "1";
格式化代码如下：
zcat Homo_sapiens.GRCh38.95.chr.gtf.gz |grep -
w havana |perl -F"\t" -
alne '{next if $F[2] ne "exon";@a=split(/\"/,$F[8]); print join("\t", $a[1],$a[5],$a[9],$F[0],$F[3],$F[4],$F[4]-
$F[3]+1) }' > g2t2exon.txt
zcat Homo_sapiens.GRCh38.95.chr.gtf.gz |grep -
w havana |perl -F"\t" -
alne '{next if $F[2] ne "CDS";@a=split(/\"/,$F[8]); print join("\t",$ a[1],$a[5],$a[9],$F[0],$F[3],$F[4],$F[4]-$F[3]+1) }' > g2t2cds.txt
结果如下：
ENSG00000223972 ENST00000456328 11118691222 7359
ENSG00000223972 ENST00000456328 21126131272 1109
ENSG00000223972 ENST00000456328 31132211440 91189
ENSG00000223972 ENST00000450305 11120101205 748
ENSG00000223972 ENST00000450305 21121791222 749
ENSG00000223972 ENST00000450305 31126131269 785
ENSG00000223972 ENST00000450305 41129751305 278
ENSG00000223972 ENST00000450305 51132211337 4154
ENSG00000223972 ENST00000450305 61134531367 0218
ENSG00000227232 ENST00000488147 11295342957 037
可以看到，一个基因会有多个转录本，然后每个转录本有多个外显子。

不同的转录本的外显子有重叠。

值得注意的是很多gtf里面记录的基因，在cDNA序列文件里面是不存在的，不过这个不影响我们的任务。

我上面两个代码分别得到的是外显子和CDS的坐标，后来发现CDS才是正确的，因为并不是所有的外显子序列都会出现在cDNA序列里面。

首先基因组坐标转为转录本坐标
接下来需要写脚本把我们转录本融合位点那个基因组坐标，转为其转录本的相对坐标，这个时候普通的shell脚本已经无能为力，需要python或者perl这样的编程语言啦，就是把我们的pos.txt文件和自己制作的g2t2exon.txt关联起来，所以需要使用关联数组这个东西。

当然，也可以使用大家最擅长的R语言展示咯。

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('input.txt')
head(a)
b=read.table('g2t2cds.txt')
colnames(b)=c('gene','transcript','exon','chr','start','end','exo n_length')
head(b)
library(stringr)
re=NULL
tmp = apply(a,1,function(x){
# x=a[1,]
g=str_split(x[1],'[.]')[[1]][1]
pos=as.numeric(str_split(x[2],':')[[1]][2])
info=b[b[,1]==g,] ## one gene might has more than 1 trans cripts.
lapply(split(info,info[,2]), function(y){
## each transcript
apply(y,1,function(z){
## each exon.
if(z[5] <= pos & z[6] >= pos ){
#print(c(z,pos))
new = sum(as.numeric(y[y[,3] < as.numeric(z[3]),7])) + pos - as.numeric(z[5])+1
print(c(z,pos,new))
re <<- rbind(re,c(z,pos,new))
}
}) ## each exon
}) ## each transcript
}) ## each position
colnames(re)=c('gene','transcript','exon','chr','start','end','ex on_length','pos','new_pos')
代码非常复杂，需要一定编程水平才能理解。

1 ENSG00000121879 ENST00000643187 2
2
3 17923409
4 17 9235098 100
5 179234680 3712
2 ENSG00000074800 ENST00000646370 1
3 1 8861007 8861 429 423 8861330 1738
3 ENSG00000074800 ENST00000647408 12 1 8861002 8861 429 428 8861330 1683
4 ENSG00000105976 ENST00000397752 1 7 116672392 11 6672577 186 116672464 73
5 ENSG0000010597
6 ENST00000456159 1
7 116672390 116 672577 18
8 116672464 75
6 ENSG00000146648 ENST00000420316
7 7 55154011 551 54152 142 55154029 1010
7 ENSG00000146648 ENST00000455089 6 7 55154011 551 54152 142 55154029 888
8 ENSG00000077782 ENST00000526570 1 8 38427921 384 30820 2900 38428753 833
9 ENSG00000037474 ENST00000502932 2 5 6603869 66042 76 408 6604266 502
如上所述，融合位点的基因组坐标，成功转为了其对应的转录本序列坐标。

比如ENST00000643187这个转录本的坐标是
ENST00000643187.1cds chromosome:GRCh38:3:179148574 :179235098:1 gene:ENSG00000121879.4 gene_biotype:protein_c oding
而融合位点在179234680，如果纯粹的使用它减去转录本起始坐标后是 86106 , 包含了大量的intron序列，所以需要找到其精准的外显子坐标。

比如这里的第22个外显子坐标是3 179234094 179235098，得到 586 的长度，再加上这个转录本前面的所有CDS的长度之和，最后是 3712 ，就是该融合位点的转录本坐标啦。

然后根据转录本坐标及转录本序列获取
这个代码也不简单，需要读取我们下载好的cDNA序列文件：
rm(list = ls())
options(stringsAsFactors = F)
load(file = 're.Rdata')
fa=readLines('Homo_sapiens.GRCh38.cds.all.fa.gz')
fa=paste0( fa ,collapse = '')
fa=strsplit(fa,'>')[[1]]
fa=fa[-1]
tid=str_split(fa,'[.]',simplify = T)[,1]
seq=str_split(fa,']',simplify = T)[,2]
head(tid)
head(seq)
re=as.data.frame(re)
re=re[re$transcript %in% tid,]
re$tseq=seq[match(re$transcript,tid)]
head(re)
re$up=unlist(lapply(1:nrow(re),function(i){
#i=1
new_pos=as.numeric(re[i,9])
l=nchar(re[i,10])
if(new_pos>500){
return(substring(re[i,10],new_pos-500,new_pos)) }else{
return(substring(re[i,10],0,new_pos))
}
}))
re$down=unlist(lapply(1:nrow(re),function(i){
#i=1
new_pos=as.numeric(re[i,9])
l=nchar(re[i,10])
if((l-new_pos) >500){
return(substring(re[i,10],new_pos,new_pos+500)) }else{
return(substring(re[i,10],new_pos,l))
}
}))
判断前面转换好的转录本坐标是否扩充500后会越界，然后选取扩充500bp的序列即可。

当然，这个结果请一定要去IGV检查哦，我这里只是提供一种思路。

值得注意的是转录本其实还有正负链信息是需要考虑的。

STAR-Fusion 提供工具组建融合后的转录本
其实并不需要自行写脚本进行探索，研究一下STAR-Fusion '--denovo_reconstruct' parameter. This requires that you include the '--FusionInspector inspect|validate' setting.
两个参数稍微有点区别
•If FusionInspector 'inspect' mode is invoked, then only the fusion-evidence reads are de novo assembled.
•If FusionInspector 'validate' mode is selected, then all reads aligned to the fusion gene contigs are assembled.
同理，上述代码可以使用perl或者python来完成，小伙伴赶快试试看吧！
历史题目：
生信编程直播第0题-生信编程很简单！
生物信息学技能面试题(第1题)-人类基因组的外显子区域到底有多长
生物信息学技能面试题(第2题)-探索人类基因组序列
生物信息学技能面试题(第3题)-探索人类基因组注释文件
生物信息学技能面试题(第4题)-多个同样的行列式文件合并起来生物信息学技能面试题(第5题)-根据GTF画基因的多个转录本结构
生物信息学技能面试题(第6题)-下载最新版的KEGG信息，并且解析好
生信编程直播第七题：写超几何分布检验！
生信编程直播第八题：ID转换大全
生信编程直播第9题-根据指定染色体及坐标得到参考碱基
生信编程直播第10题：根据指定染色体及坐标得到位置信息
生信编程直播第11题：把文件内容按照染色体分开写出。