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

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

首先需要理解:基因,转录本(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 SpanningFragCount SpliceType LeftGene LeftBreakpoint RightGene RightBreakpoint LargeAnchorSupport FFPM LeftBreakDinuc LeftBreakEntropy RightBreakDinuc RightBreakEntropy annots
THRA--AC090627.1 27 93 ONLY_REF_SPLICE THRA^ENSG00000126351.8 chr17:38243106:+ AC090627.1^ENSG00000235300.3 chr17:46371709:+ YES_LDAS 23875.8456 GT 1.8892 AG 1.9656 ["CCLE","FA_CancerSupp","INTRACHROMOSOMAL[chr17:8.12Mb]"]

可以看到列比较多,值得我们关心的是转录本融合的左右两个转录本的ID及其融合位点在基因组的坐标,如下:

THRA^ENSG00000126351.8 chr17:38243106:+ 
AC090627.1^ENSG00000235300.3 chr17:46371709:+

上面的坐标可以无限多,文件命名为 pos.txt 吧,后面写脚本需要用到,值得注意的是作者软件举例应该是hg19的基因组坐标,目前几乎都是hg38了,所以后面我脚本也是基于hg38的。

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

我们很容易拿到各个转录本的基因组坐标,但是融合点的基因组坐标不能简单对应到转录本cDNA序列里面坐标。我们的突破点,就是找到融合点的基因组坐标到底对应到转录本cDNA序列的哪个位置。

物种的基因及其注释文件的下载

这里首选:https://asia.ensembl.org/info/data/ftp/index.html

wget ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/cds/Homo_sapiens.GRCh38.cds.all.fa.gz
wget ftp://ftp.ensembl.org/pub/release-95/gtf/homo_sapiens/Homo_sapiens.GRCh38.95.chr.gtf.gz

首先解析基因组注释文件,就是gtf的:

 1 havana exon 11869 12227 . + . gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; 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 1 1 11869 12227 359
ENSG00000223972 ENST00000456328 2 1 12613 12721 109
ENSG00000223972 ENST00000456328 3 1 13221 14409 1189
ENSG00000223972 ENST00000450305 1 1 12010 12057 48
ENSG00000223972 ENST00000450305 2 1 12179 12227 49
ENSG00000223972 ENST00000450305 3 1 12613 12697 85
ENSG00000223972 ENST00000450305 4 1 12975 13052 78
ENSG00000223972 ENST00000450305 5 1 13221 13374 154
ENSG00000223972 ENST00000450305 6 1 13453 13670 218
ENSG00000227232 ENST00000488147 1 1 29534 29570 37

可以看到,一个基因会有多个转录本,然后每个转录本有多个外显子。不同的转录本的外显子有重叠。值得注意的是很多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','exon_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 transcripts.
 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','exon_length','pos','new_pos')

代码非常复杂,需要一定编程水平才能理解。

1 ENSG00000121879 ENST00000643187 22 3 179234094 179235098 1005 179234680 3712
2 ENSG00000074800 ENST00000646370 13 1 8861007 8861429 423 8861330 1738
3 ENSG00000074800 ENST00000647408 12 1 8861002 8861429 428 8861330 1683
4 ENSG00000105976 ENST00000397752 1 7 116672392 116672577 186 116672464 73
5 ENSG00000105976 ENST00000456159 1 7 116672390 116672577 188 116672464 75
6 ENSG00000146648 ENST00000420316 7 7 55154011 55154152 142 55154029 1010
7 ENSG00000146648 ENST00000455089 6 7 55154011 55154152 142 55154029 888
8 ENSG00000077782 ENST00000526570 1 8 38427921 38430820 2900 38428753 833
9 ENSG00000037474 ENST00000502932 2 5 6603869 6604276 408 6604266 502

如上所述,融合位点的基因组坐标,成功转为了其对应的转录本序列坐标。

比如 ENST00000643187 这个转录本的坐标是

ENST00000643187.1 cds chromosome:GRCh38:3:179148574:179235098:1 gene:ENSG00000121879.4 gene_biotype:protein_coding

而融合位点在 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的序列即可。

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

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.

#

Comments are closed.