本文主要介绍TEtranscripts，常用于转座子研究。
官网：https://github.com/mhammell-laboratory/TEtranscripts/
Jin Y, Tam OH, Paniagua E, Hammell M. (2015) TEtranscripts: a package for including transposable elements in differential expression analysis of RNA-seq datasets. Bioinformatics. 31(22):3593-9.
简介
转座子其实属于重复序列中的一个特殊类别，由于重复序列的研究起来比较复杂，所以一直以来研究这方面的人并不多。今天我们就来说一说如何分析差异转座子，其实从本质上看与分析差异基因表达类似，把转座子当成基因，然后统计reads，最后利用差异软件如DESeq2做差异分析。但是由于转座子具有重复序列，所以序列间相似性很高，所以比对到该区域的reads有很多是mutil-alignments，也就是说同一条read会有好多个比对到的地方。如果简单的统计reads，很多 read被重复计数，如果只用unique-alignments的reads，又会丢失很多信息。所以研究人员开发了相应的分析软件 — TEtranscripts。该软件在定量时，考虑了比对到转座子区域的reads数目(给予每一条mutil-alignments的reads一个权重)，转座子的长度，fragment长度等信息，使用Expectation maximization方法确定reads具体属于哪一个转座子。
从流程图可以看出该软需要三种输入文件，样本bam、基因组gtf、转座子gtf，包含两个子命令TEtranscripts、TEcount，前一个命令会输出表达值和差异结果，后一个只是统计表达值。统计reads的时候，软件有两种模式可选分别为multi、uniq，默认是multi。如果选择uniq模式，mutil-alignments的reads就会被丢弃不用于计数。一般选择默认情况更好。该软件由python编写，安装和使用都很方便。
安装
方法一：克隆
cd ~/software/
git clone https://github.com/mhammell-laboratory/TEtranscripts
cd ~/software/TEtranscripts
python setup.py install --user
TEtranscripts --version
#输出：TEtranscripts 2.2.3
方法二：压缩包
cd ~/software/
wget https://files.pythonhosted.org/packages/a1/e7/9f6a4d3d466374c22fe732d79a6ae40fe86384182e97085b2b9a343b00db/TEtranscripts-2.2.3.tar.gz
tar -xzvf TEtranscripts-2.2.3.tar.gz
# 建议把解压后的文件夹名称修改为TEtranscripts
cd ~/software/TEtranscripts/
python3 setup.py install --user
TEtranscripts --version
TEtranscripts -h
PS：
- 请先去官网学习了解基本信息
- 查看一下安装的setup.py中是version='2.2.3'，老版本使用python3会有很多error和bug
GTF下载
这一块一定要确保两个GTF文件是同一种类型的，可以查看前几行gene_id是不是同种类型的，不要将genecode和ensenbal的混用，容易导致无法定量
### 设置保存GTF文件的目录
cd ~/software/reference/GRCh38/
#### Gene 文件
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/genes/hg38.ncbiRefSeq.gtf.gz
gunzip hg38.ncbiRefSeq.gtf.gz
#### TE文件
wget https://labshare.cshl.edu/shares/mhammelllab/www-data/TEtranscripts/TE_GTF/GRCh38_GENCODE_rmsk_TE.gtf.gz
# 解压
gunzip *.gz
# 重命名
mv hg38.ncbiRefSeq.gtf GRCh38_genecode.gtf
该网站还可以下载很多其他物种的TF文件（https://labshare.cshl.edu/shares/mhammelllab/www-data/TEtranscripts/TE_GTF/）
配置
由于差异分析需要用到R和DEGSeq2，如果没有安装请先安装这两个软件
# 安装R和biocmanager（用于管理R包）
conda install -c conda-forge r-base
conda install -c bioconda bioconductor-biocmanager
# 启动R
R
# 使用BiocManager安装 DESeq2
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")
library("DESeq2")
#输出一大段且不报错说明安装成功
# 退出R环境
q()
参数
具体且详细的参数设置请去看官网！！！
下面对TEtranscripts命令的常用且重要的几个参数进行说明（TEcount类似，请看官网）
TEtranscripts
-t treatment sample [实验组样本，可以有多个重复]
-c control sample [对照组样本，可以有多个重复]
--GTF [基因组注释参考文件]
--TE [转座子注释参考文件]
--format [输入文件类型：bam或sam，默认是bam]
--stranded [是否为链特异性以及链特异性的方向]
no - 非链特异性文库，默认
forward - "Second-strand" cDNA 文库，来源于第二条反向链(例：QIAseq stranded)
reverse - "First-strand" cDNA 文库，来源于第一条正向链(例：Illumina TruSeq stranded)
--sortByPos [输入文件是按照染色体位置排序后的]（如果不是则不用加）.
--project [输出文件的前缀名称]
--outdir [输出文件夹路径]
--mode [TE计数模式]
uniq (唯一比对)
multi (在所有队列中分配，默认且推荐用multi).
-p | --padj [FDR值，默认为0.05]
-f | --foldchange [变化倍数默认为1]
--DESeq [使用DESeq做差异分析，默认是DESeq2]
-n | --norm [标准化方式]
DESeq_default (默认),
TC (total annotated read counts),
quant (quantile normalization).
使用
我的bam文件是bulk RNA的数据经过STAR比对得到后的文件，开发者建议用STAR进行比对，因为这个包是专门为RNA-seq开发的。（其他的对比软件得到的结果也可以用）
### TEtranscripts差异分析
cd ~/project/TE/test/DETE
TEtranscripts --mode multi \
--format BAM \
-t ../bam/TSC_KO_rep1.bam ../bam/TSC_KO_rep2.bam \
-c ../bam/TSC_WT_rep1.bam ../bam/TSC_WT_rep2.bam \
--GTF ~/software/reference/GRCh38/GRCh38_genecode.gtf \
--TE ~/software/reference/GRCh38/GRCh38_GENCODE_rmsk_TE.gtf \
--foldchange 1 --sortByPos \
--project TSC
### TEcount 表达分析
cd ~/project/TE/test/expr
for i in TSC_KO_rep1 TSC_KO_rep2 TSC_WT_rep1 TSC_WT_rep2; do
TEcount --format BAM --mode multi -b ../bam/${i}.bam \
--GTF ~/software/reference/GRCh38/GRCh38_genecode.gtf \
--TE ~/software/reference/GRCh38/GRCh38_GENCODE_rmsk_TE.gtf \
--sortByPos --project ${i}
done
结果
DETE
├── TSC.cntTable # count矩阵
├── TSC_DESeq2.R # 差异R脚本
├── TSC_gene_TE_analysis.txt # 总的差异结果
└── TSC_sigdiff_gene_TE.txt # 显著性的差异结果
expr
├── TSC_WT_rep1.cntTable # 单样本中GENE_TE表达
├── TSC_WT_rep2.cntTable
├── TSC_KO_rep1.cntTable
└── TSC_KO_rep2.cntTable
-
全部基因的表达和差异变化信息
TSC.cntTable
- 卡阈值后的差异基因信息（可以看到我的数据有7461个显著差异的）
后续可以对差异的TE进行功能富集等分析
-
TEcounts得到的每个文件都包含两列：一列名称一列表达
TSC_KO_rep1.cntTable
注意事项：
- 不要残留老版本的TEtranscripts，卸载干净再从新安装，或者直接删除所有相关的包
- 两个gtf文件的染色体名称要一致，例如都加chr，或者都不加
- TF的GTF文件的第四五列要是数值型的，不能是科学计数法
- 环境冲突/版本不兼容，如果各种文件都确保了没有问题但仍有奇怪的报错，建议创建虚拟环境然后运行
# 指定python版本创建TE_Transcripts虚拟环境
cd ~/software/tools
conda create -n TE_Transcripts python=3.10
# 激活
conda activate TE_Transcripts
# 安装 TEtranscripts
git clone https://github.com/mhammell-laboratory/TEtranscripts
cd ./TEtranscripts
python setup.py install --user
# 检验安装是否成功
~/software/tools/TEtranscripts/bin/TEtranscripts --version
~/software/tools/TEtranscripts/bin/TEtranscripts -h
# 安装R和DESeq2，进行后续分析....
- 官网issue部分可以搜索已知的报错，如果还是没有解决可以直接发布新issue问作者，回复超快！👍
欢迎大家评论交流！
（每帖分享：日升日落，总有黎明！）
