写在前面,最近为了解决这个问题,看了很多软件,要么是单细胞技术,要么是物种限制,总之花了挺多时间。本来打算使用SQuIRE,但是卡在了BUILD_refGene.txt(基因注释和TE文件的格式上),只得作罢,改用这个TEtranscripts.
- 下载
wget https://github.com/mhammell-laboratory/TEtranscripts/archive/refs/tags/2.2.3.zip
解压并进入文件夹
tar -zxvf TEtranscripts-2.2.3.tar.gz && cd TEtranscripts-2.2.3
安装
python setup.py install --user
- 基因注释文件和TE注释文件
在github的issue中找到了作者提供的脚本,详情查看该链接
https://github.com/mhammell-laboratory/TEtranscripts/issues/113
下载该脚本(本软件需要特殊格式的TE注释文件)
wget https://labshare.cshl.edu/shares/mhammelllab/www-data/TEtranscripts/TE_GTF/makeTEgtf.pl.gz
gunzip makeTEgtf.pl.gz
利用该脚本产生合适的gtf的TE注释文件(根据自己文件指定相应列数)
./makeTEgtf.pl -c 2 -s 3 -e 4 -o 5 -t 1 -f 6 Maize_TE.SAF > Maize_TE.gtf
SAF格式如下
SAF.png
3.使用star进行比对
3.1下载并且安装
wget https://github.com/alexdobin/STAR/archive/2.7.10b.tar.gz
tar -xzf 2.7.10b.tar.gz
cd STAR-2.7.10b
cd ./source
make STAR
3.2 建立索引
STAR --runThreadN 6 --runMode genomeGenerate --genomeDir ./ --genomeFastaFiles ref.fasta --sjdbGTFfile ref.gtf --sjdbOverhang 99
参数解释
sjdbOverhang: 在多数情况下,默认值 100 的作用与理想值类似。
3.3 进行比对
STAR --genomeDir Path/to/genome_index --readFilesIn /Path/to/RNA.fq.gz --outFileNamePrefix rep1 --sjdbScore 2 --runThreadN 8 --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --winAnchorMultimapNmax 100 --outFilterMultimapNmax 100
参数解释
sjdbScore: extra alignment score for alignmets that cross database junctions
winAnchorMultimapNmax:max number of loci anchors are allowed to map to
outFilterMultimapNmax:maximum number of loci the read is allowed to map to
4.利用 TEtranscripts 进行定量和差异分析, 可以直接得到DESeq2的结果
TEtranscripts -t treatment.bam -c control.bam --GTF SollycM82_genes_v1.1.0_1.gtf --TE M82_TE.gtf --mode multi --project name
4.1 针对每一个样本进行定量,然后就和基因的DESeq2分析一样了
TEcount --sortByPos --format BAM --mode multi -b sample1.sortedByCoord.out.bam --GTF maize_genes.gtf --TE Maize_TE.gtf --project sample1