本文主要介绍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具体属于哪一个转座子。

TEtranscripts flow chart

从流程图可以看出该软需要三种输入文件，样本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：

1. 请先去官网学习了解基本信息
2. 查看一下安装的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/）

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

1. 全部基因的表达和差异变化信息

   

   TSC.cntTable

TSC_gene_TE_analysis.txt

2. 卡阈值后的差异基因信息（可以看到我的数据有7461个显著差异的）

TSC_sigdiff_gene_TE.txt

后续可以对差异的TE进行功能富集等分析

3. 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问作者，回复超快！👍

欢迎大家评论交流！
（每帖分享：日升日落，总有黎明！）