自制转座子GTF

转座子数据库

  上一篇《TEtranscripts：转座子元件差异表达》介绍了如何用TEtranscripts软件来做转座子差异分析。但是，在使用软件前我们需要准备好必需的两个GTF文件，基因的GTF比较常见好多网站如gencode都可以下载到，这里就不再赘述了。
  哪里能下载到转座子GTF？可能没有很直接的答案，至少本人是没有找到可以直接下载到转座子GTF的数据库。本来把希望都寄托在Repbase数据库上，作为最常用的DNA重复序列数据库，里面也是包含了比较全面的转座子序列。很遗憾😞😞😞，这个数据库现在已经收费了。既然没有直接方法，只能绕一个弯，找一个不收费的数据库下载数据，自制一个转座子的GTF。

下载数据

  UCSC网站大家应该不陌生，这里就不介绍了，庆幸的是该网站收录了来自Repbase和RepeatMasker两个数据库的重复序列，由于Repbase从2010年开始收费，所以UCSC对repeatmasking注释的最后一次更新也停在了2010年。既然我们想要免费的数据，那就不能有什么要求了，毕竟现在没有更好的白漂方式了。

  从上面可以知道，UCSC收录的重复序列的种类，不知道这些种类里面哪些是转座子的可以戳这里TEtranscripts：转座子元件差异表达。

# 下载
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/database/rmsk.txt.gz
gzip -cd rmsk.txt.gz >hg38_rmsk.txt

  我这里下载的是hg38的重复序列，如果你需要下载小鼠的重复序列可以将上面的网址里面的hg38替换为mm10即可。当然，还有hg19、mm9也是可以的。下载的文件内容如下：

制作GTF

  可以看到文件里面包含了我们需要的信息，重复序列的坐标、名称、家族名、类名，根据这些信息我们就可以制作成GTF文件了，下面我们就可以开干了：

awk '{OFS=FS="\t";a[$11]++;print$6,"hg38_rmsk","exon",$7,$8,$9,$10,".",$11,$11"_dup"a[$11]-1,$13,$12}' hg38_rmsk.txt | \
awk '{OFS=FS="\t";gsub("_dup0","",$10);print$1,$2,$3,$4,$5,$6,$7,$8,"gene_id ""\""$9"\"; ""transcript_id ""\""$10"\"; ""family_id ""\""$11"\"; ""class_id ""\""$12"\";"}' >hg38_rmsk.gtf

结果如下：

  生成gtf还是很容易的，我这里就直接使用了awk命令来完成了，对于这样简单的需求使用awk很方便就可以实现。作为linux三剑客的awk功能可是相当强大的，用好了可以解决不少问题，本人列举了一些使用场景，感兴趣的可以戳这里Linux： 数据处理三剑客之awk，文件去重：awk的高级用法。
  重点来了！可能有人会疑惑，现在生成的GTF里面不全是转座子呀，是不是应该去除其他不是转座子的重复序列呢？不应该删除，GTF应该尽量全面，这是为了防止发生类似文库补偿的情况。我们都知道重复序列本身就有重复，重复序列之间的相似性也比较高，例如一条read原本回帖到A、B两个地方，如果我们准备的GTF里面恰巧没有A的信息，那这条read就会被直接计数给B了，这样就会造成结果的偏向性。所以，让GTF信息全面可以提高结果的准备性，因为TEtranscripts软件在计数时，首先会根据算法判断reads(同时回帖到多个地方)最有可能属于哪个地方。

结束语

  至此，我们分析转座子的差异应该没啥问题了。当然，分享的内容也都是本人自己学习后的认知，或许会有哪些说的不对的地方，希望知道的读者能够予以指出。。。