TransposonPSI——转座子分析的入门自学
最近需要做转座子分析，查找发现可以使用 TransposonPSI 来进行分析。但是登陆官网，该软件 update 时间为 2013 年，但是因为时间紧迫，暂时还没有进行其他方法的调研，所以先选用该软件进行了分析。
一、TransposonPSI 安装及使用
1. TransposonPSI 安装
官网： http://transposonpsi.sourceforge.net
下载地址：https://sourceforge.net/projects/transposonpsi/
压缩包非常小，只有 10M 左右，解压后修改主角本 transposonPSI.pl 中三个软件的路径（blastall, formatdb, blastpgp），即可食用。
目录结构：
README docs/ PerlLib/ scripts/ transposon_ORF_lib/ transposon_PSI_LIB/ misc/ transposonPSIcreate/ TransposonWeb/ transposonPSI.pl
test/
2. TransposonPSI 使用入门
直接进入 test 目录，执行 runMe.sh 即可进行测试，非常简单。查看 runMe.sh 发现，输入文件是我们需要进行分析的数据序列，nuc 表示核酸序列，prot 表示蛋白序列。
if [ -e target_test_genome_seq.fasta.gz ] && ! [ -e target_test_genome_seq.fasta ] then gunzip target_test_genome_seq.fasta.gz fi ../transposonPSI.pl target_test_genome_seq.fasta nuc
二、TransposonPSI 流程解读
对 transposonPSI.pl 进行 Linux 脚本复现
cd /Transposon/div_step/ if [ -d tmp ] then rm -rf tmp fi mkdir tmp cd tmp ln -s ../target_test_genome_seq.fasta /software/blast-2.2.26/bin/formatdb -i target_test_genome_seq.fasta -p F cd /Transposon/div_step/tmp TPSI_list=( cacta DDE_1 gypsy hAT helitron ISa ISb isc1316 line ltr_Roo mariner_ant1 mariner MuDR P_element piggybac TY1_Copia TyrRecombinaseCrypton ) for int in {0..16} do name=${TPSI_list[$int]} /software/blast-2.2.26/bin/blastall -i /Transposon/TransposonPSI_08222010/transposon_PSI_LIB/$name.refSeq -d target_test_genome_seq.fasta -p psitblastn -R /Transposon/TransposonPSI_08222010/transposon_PSI_LIB/$name.chk -F F -M BLOSUM62 -t -1 -e 1e-5 -v 10000 -b 10000 >target_test_genome_seq.$name.psitblastn /Transposon/TransposonPSI_08222010/scripts/BPbtab </Transposon/div_step/tmp/target_test_genome_seq.$name.psitblastn> /Transposon/div_step/tmp/target_test_genome_seq.$name.psitblastn.btab done cat /Transposon/div_step/tmp/*btab | sort -rn -k13 >/Transposon/div_step/target_test_genome_seq.TPSI.allHits cd /Transposon/div_step/ perl /Transposon/TransposonPSI_08222010/scripts/TBLASTN_hit_chainer.pl target_test_genome_seq.TPSI.allHits btab >target_test_genome_seq.TPSI.allHits.chains perl /Transposon/TransposonPSI_08222010/scripts/TPSI_btab_to_gff3.pl target_test_genome_seq.TPSI.allHits.chains >target_test_genome_seq.TPSI.allHits.chains.gff3 perl /Transposon/TransposonPSI_08222010/scripts/TBLASTN_hit_chainer_nonoverlapping_genome_DP_extraction.pl target_test_genome_seq.TPSI.allHits.chains >target_test_genome_seq.TPSI.allHits.chains.bestPerLocus perl /Transposon/TransposonPSI_08222010/scripts/TPSI_chains_to_gff3.pl target_test_genome_seq.TPSI.allHits.chains.bestPerLocus >target_test_genome_seq.TPSI.allHits.chains.bestPerLocus.gff3
1. 格式化序列数据库
这是 blast 比对的首要步骤，这里我就不多介绍了，详细的参数和使用说明很多大佬都有介绍，使用时百度即可。
/software/blast-2.2.26/bin/formatdb -i target_test_genome_seq.fasta -p F
2. 数据库列表准备
TPSI_list=(
cacta
DDE_1
gypsy
hAT
helitron
ISa
ISb
isc1316
line
ltr_Roo
mariner_ant1
mariner
MuDR
P_element
piggybac
TY1_Copia
TyrRecombinaseCrypton
)
以上列表为各类转座子名称，它们存在于 transposon_PSI_LIB/ 目录中，每一种数据库有三种格式：refSeq，chk，chkp
3. 序列与各数据库进行比对
/software/blast-2.2.26/bin/blastall \
-i /Transposon/TransposonPSI_08222010/transposon_PSI_LIB/$name.refSeq \
-d target_test_genome_seq.fasta \
-p psitblastn -R /Transposon/TransposonPSI_08222010/transposon_PSI_LIB/$name.chk \
-F F \
-M BLOSUM62 \
-t -1 \
-e 1e-5 \
-v 10000 \
-b 10000 \
>target_test_genome_seq.$name.psitblastn
特殊参数
-R PSI-TBLASTN checkpoint file [File In] Optional
得到比对结果。
4. BPbtab
/Transposon/TransposonPSI_08222010/scripts/BPbtab \
</Transposon/div_step/tmp/target_test_genome_seq.$name.psitblastn> \
/Transposon/div_step/tmp/target_test_genome_seq.$name.psitblastn.btab
将比对结果转化为 btab 格式：
1 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 752 784 637 735 81.8 81.8 130 54.5 0 Plus 117619 1e-09 2 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 751 791 769 891 68.3 78.0 126 52.9 0 Plus 117619 4e-09 3 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 753 781 37440 37526 89.7 89.7 124 52.2 0 Plus 117619 7e-09 4 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 752 784 622 720 75.8 75.8 118 49.8 0 Plus 117619 3e-08 5 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 753 790 658 771 68.4 71.1 118 49.8 0 Plus 117619 3e-08 6 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 752 845 41583 41933 29.7 39.8 117 49.5 0 Plus 117619 4e-08 7 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 753 789 664 774 70.3 73.0 116 49.1 0 Plus 117619 6e-08 8 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 752 785 649 750 76.5 79.4 115 48.7 0 Plus 117619 7e-08 9 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 750 838 37422 37697 33.3 46.5 115 48.7 0 Plus 117619 7e-08 10 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 756 800 775 909 55.6 64.4 114 48.3 0 Plus 117619 1e-07 11 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 751 784 625 726 73.5 73.5 111 47.2 0 Plus 117619 2e-07 12 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 752 789 715 873 58.5 60.4 111 47.2 0 Plus 117619 2e-07 13 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 753 784 682 810 62.8 62.8 108 46.0 0 Plus 117619 5e-07 14 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 753 784 706 831 61.9 61.9 103 44.1 0 Plus 117619 2e-06 15 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 751 784 577 696 62.5 62.5 102 43.7 0 Plus 117619 3e-06
BPtab 是一个Blast输出解析器， 脚本将 WU-BLAST 或 NCBI-BLAST 输出文件解析为BTAB格式，其中每个 HSP 报告为带有制表符分隔字段的单行。
5. 统计比对结果
cat /Transposon/div_step/tmp/*btab | sort -rn -k13 >/Transposon/div_step/target_test_genome_seq.TPSI.allHits
1 TY1_Copia 1643 PSITBLASTN target_test_genome_seq.fasta genome 511 1530 4007 7159 27.2 46.4 1826 707 0 Plus 117619 0 2 helitronORF 1321 PSITBLASTN target_test_genome_seq.fasta genome 663 1175 7497 9239 53.4 64.2 1633 633 0 Plus 117619 0 3 Crypton 457 PSITBLASTN target_test_genome_seq.fasta genome 1 456 85676 87118 40.7 57.7 1148 446 0 Plus 1.18E-139 4 TY1_Copia 1643 PSITBLASTN target_test_genome_seq.fasta genome 1149 1641 52051 53538 35.3 56.2 1138 442 0 Plus 117619 1.00E-129 5 helitronORF 1321 PSITBLASTN target_test_genome_seq.fasta genome 998 1321 35448 36542 56 66 1086 422 0 Plus 117619 1.00E-125 6 gypsy 1463 PSITBLASTN target_test_genome_seq.fasta genome 574 1018 111538 112860 26.7 46.7 1012 394 0 Plus 1.18E-109 7 cacta 1264 PSITBLASTN target_test_genome_seq.fasta genome 752 784 637 735 81.8 81.8 130 54.5 0 Plus 1.18E-04
BTAB 格式的具体内容并为完全掌握，暂时不提。
6. 共线性 HSPs 关联
perl /Transposon/TransposonPSI_08222010/scripts/TBLASTN_hit_chainer.pl \
target_test_genome_seq.TPSI.allHits btab \
>target_test_genome_seq.TPSI.allHits.chains
collinear HSPs are chained together into larger chains (more complete element regions).
将共线性的 HSP 连接在一起，形成 larger chains，例如下面的文件，会将线性相关的放在一起
1 #Chain Crypton 167-308 genome 23349-23753 + 46.3 2 Crypton 457 PSITBLASTN target_test_genome_seq.fasta genome 167 308 23349 23753 32.6 47.9 210 85.5 3 4 #Chain Crypton 1-257 genome 37524-38356 + 92.7 5 Crypton 457 PSITBLASTN target_test_genome_seq.fasta genome 1 73 37524 37742 42.5 57.5 160 66.2 6 Crypton 457 PSITBLASTN target_test_genome_seq.fasta genome 74 257 37742 38356 33.5 50.0 297 118 7 8 #Chain Crypton 52-456 genome 40190-41483 + 148.0 9 Crypton 457 PSITBLASTN target_test_genome_seq.fasta genome 52 189 40190 40600 37.9 55.7 258 103 10 Crypton 457 PSITBLASTN target_test_genome_seq.fasta genome 182 456 40641 41483 34.8 48.2 401 159
7. 将 larger chains 转 gff3 格式
perl /Transposon/TransposonPSI_08222010/scripts/TPSI_btab_to_gff3.pl \
target_test_genome_seq.TPSI.allHits.chains \
>target_test_genome_seq.TPSI.allHits.chains.gff3
8. best chains ：提取分数最高的 chains
perl /Transposon/TransposonPSI_08222010/scripts/TBLASTN_hit_chainer_nonoverlapping_genome_DP_extraction.pl \
target_test_genome_seq.TPSI.allHits.chains \
>target_test_genome_seq.TPSI.allHits.chains.bestPerLocus
9. best chains 转 gff3 格式
perl /Transposon/TransposonPSI_08222010/scripts/TPSI_chains_to_gff3.pl \
target_test_genome_seq.TPSI.allHits.chains.bestPerLocus \
>target_test_genome_seq.TPSI.allHits.chains.bestPerLocus.gff3
即为最终结果。
