复现Tripal Synview功能模块

By masiyi, 11 May, 2024

Forums

Researches

本模块依然在91的师姐的rhoddb-test容器中进行测试

一、模块的下载和安装

在GitHub上面通过链接下载或者下载压缩包到本地（https://github.com/tripal/tripal_synview/blob/master/README.md）

二、所需数据收集和前期处理

1.数据收集

这里测试以RPGD上的马樱杜鹃V2版本的基因组fasta、注释gff3和蛋白序列fasta，圆叶杜鹃基因组fasta、注释gff和蛋白序列fasta为测试数据。均可在RPGD上下载。

2.数据处理

由于此模块是将最终计算的种内（种间）基因共线性区块结果给可视化了，所以还需要得到这个共线性结果。在这里我们使用MCScanX，也是从GitHub上可以下载。

2.1gff文件处理成所需格式

MCScanX需要gff文件和blast结果文件作为输入文件，由于它所需要的gff文件和标准gff文件不一样，所以我们需要先对gff文件进行处理（这里以马樱杜鹃的gff为例），提取它需要的列生成一个新的gff文件，在这里我们只要cds，首先先提取所有cds，命令如下：

grep -w "CDS" Rhododendron_delavayi.gene.gff3 > cds.gff3

接着提取我们需要的四列并且写为一个新的gff文件，这里用了一个awk命令的脚本（因为这个awk命令有点长），脚本 tiquGFF.awk 如下（这个是否选择冒号要根据具体文件进行调整，记住第九列我们只要基因id就可以了）：

{  
    # 提取第九列的内容  
    attr = $9;  
      
    # 查找ID=并提取其后的值直到下一个冒号  
    start = index(attr, "ID=");  
    if (start != 0) {  
        end = index(substr(attr, start+3), ":");  
        if (end == 0) {  
            # 如果找不到冒号，则取到字符串末尾  
            id = substr(attr, start+3);  
        } else {  
            # 否则取到冒号之前  
            id = substr(attr, start+3, end-1);  
        }  
    } else {  
        # 如果没有ID字段，id变量设置为空字符串  
        id = "";  
    }  
      
    # 打印第一、四、五列和提取的ID值，以制表符分隔  
    print $1, id, $4, $5;  
}

命令行输入 `awk -f tiquGFF.awk cds.gff3 > Rhododendron_delavayi_MCScanX_v1.gff3` 即可调用脚本执行并且输出我们需要的最终版本的 Rhododendron_delavayi_MCScanX_v1.gff3 文件，此文件大致格式如下，其中第一列是染色体序号、第二列是基因id、三四列是起止位置：

rhdel01 Rhdel01G0000100.1 31115 31164
rhdel01 Rhdel01G0000100.1 31250 31325
rhdel01 Rhdel01G0000100.1 33214 33282
rhdel01 Rhdel01G0000100.1 33411 33576

以同样的方法我们可以得到圆叶杜鹃的 Rhododendron_williamsianum_MCScanX_v1.gff 文件

2.2blast种内种间关系

这里首先下载并安装ncbi的blast+，同样在GitHub上下载，在命令行输入命令 `wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.9.0/ncbi-blast-2.9.0+-x64-linux.tar.gz ` 即可，解压 `tar -xzvf ncbi-blast-2.9.0+-x64-linux.tar.gz `

接着进入到它的bin目录下可以看到它所有的命令

2.2.1建blast库

这里我们建立两个蛋白库，马樱和圆叶的，还是以马樱为例，在bin目录下输入命令

 makeblastdb -in /var/www/html/sites/default/files/raw/Rhododendron_delavayi/Rhododendron_delavayi.pep.fa -dbtype prot -out Rhododendron_delavayi.db

成功建库后我们会得到一系列索引文件，为了容易辨认和找到，我新建一个Rhododendron_delavayi_db的目录将它们放在了上级目录下面。以同样的方法构建了圆叶杜鹃的蛋白库

2.2.2比对

这里我们用圆叶杜鹃去比对马樱的蛋白库，然后再用马樱的去比圆叶的蛋白库。以圆叶杜鹃比马樱为例：

./blastp -query /var/www/html/sites/default/files/raw/Rhododendron_williamsianum/R.will.pep.fa -db /var/www/html/sites/default/files/raw/Rhododendron_delavayi_db/Rhododendron_delavayi.db -outfmt 6 -out ./rhwil_rhdel_final.blast

这里同样方法会得到两个blast结果文件，

./blastp -query /var/www/html/sites/default/files/raw/Rhododendron_delavayi/Rhododendron_delavayi.pep.fa -db /var/www/html/sites/default/files/raw/Rhododendron_williamsianum_db/Rhododendron_williamsianum.db -outfmt 6 -out ./rhdel_rhwil_final.blast

最好还是.blast结尾命名吧，方便后面使用MCScanX工具

2.2.3使用MCScanX

首先合并上面得到的gff和blast

cat Rhododendron_delavayi_MCScanX_v1.gff3 Rhododendron_williamsianum_MCScanX_v1.gff > rhdel_rhwil.gff
cat rhdel_rhwil_final.blast rhwil_rhdel_final.blast >rhdel_rhwil.blast

将这两个文件放到MCScanX下面的同一个目录下（都在 /var/www/html/sites/default/files/raw/MCScanX-master 这个里面），运行命令 MCScanX rhdel_rhwil ，还没试，感觉会报错，感觉文件名可能需要和gff或者blast结果里面的基因id相对应，完了先试试不对应行不行，不行的话再对应，开始尝试

PS:如果你在输出信息中发现 0 matches imported (xxxxx discarded), 那么一定是你的GFF文件里的基因名和blast结果的基因名不对应导致

现在运行后有问题了：

[root@f6f5e8d1a851 MCScanX-master]# ./MCScanX rhdel_rhwil
Reading BLAST file and pre-processing
Generating BLAST list
0 matches imported (1 discarded)
0 pairwise comparisons
0 alignments generated
Pairwise collinear blocks written to rhdel_rhwil.collinearity [0.368 seconds elapsed]
Writing multiple syntenic blocks to HTML files
QEFC01011942.1 cds-KAE9444543.1 834 935.html
Done! [0.030 seconds elapsed]

接着我拿这个软件中示例文件跑了一下，结果出了：

[root@f6f5e8d1a851 MCScanX-master]# ./MCScanX data/example
Reading BLAST file and pre-processing
Generating BLAST list
355520 matches imported (237427 discarded)
275 pairwise comparisons
1358 alignments generated
Pairwise collinear blocks written to data/example.collinearity [8.191 seconds elapsed]
Tandem pairs written to data/example.tandem
Writing multiple syntenic blocks to HTML files
os1.html
os10.html
os11.html
os12.html
os2.html
os3.html
os4.html
os5.html
os6.html
os7.html
os8.html
os9.html
sb1.html
sb10.html
sb2.html
sb3.html
sb4.html
sb5.html
sb6.html
sb7.html
sb8.html
sb9.html
zm1.html
zm10.html
zm2.html
zm3.html
zm4.html
zm5.html
zm6.html
zm7.html
zm8.html
zm9.html
Done! [0.730 seconds elapsed]

到这里中止了，应该是圆叶杜鹃数据不好导致的

三、使用示例文件进行流程

接着往下做就是使用MCScanX软件下面的示例数据at_vv这组数据（也记得通过上面的步骤得到.collinearity文件）

在这之前写一个at_vv.chr文件，两列，第一列是原本的染色体名称，第二列是新的染色体名称，里面要写全at和vv这两个文件的染色体，如下：

chr1    at1
chr2    at2
chr3    at3
chr4    at4
chr5    at5
chr1    vv1
chr2    vv2
chr3    vv3
chr4    vv4
chr5    vv5
chr6    vv6
chr7    vv7
chr8    vv8
chr9    vv9
chr10   vv10
chr11   vv11
chr12   vv12
chr13   vv13
chr14   vv14
chr15   vv15
chr16   vv16
chr17   vv17
chr18   vv18
chr19   vv19
chr20   vv20
chr21   vv21
chr22   vv22
chr23   vv23
chr24   vv24
chr25   vv25
chr26   vv26
chr27   vv27
chr28   vv28
chr29   vv29

现在根据要求跑一个准备文件 .block文件：

perl ./syntenyTool.pl -t mcscanx_block -c /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.chr /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.collinearity /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.gff > /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block

四、tripal所需文件的生成

下一步要查找这两个物种的organism_id，我就给它算到1和2里面吧，反正测试，用下列命令生成可以加载到chado里的文件：

[root@f6f5e8d1a851 tripal_synview-master]# perl ./syntenyTool.pl -t mcscanx_tripal -a at -b vv -c 1 -d 2 /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.chr /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.gff /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.collinearity /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block > /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block.tripal.txt

截至到目前为止，一直会出现报错：

[root@f6f5e8d1a851 tripal_synview-master]# perl ./syntenyTool.pl -t mcscanx_tripal -a at -b vv -c 1 -d 2 /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.chr /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.gff /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.collinearity /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block > /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block.tripal.txt
[ERR]undef ref2 chr vv at1

应该还是输入文件不太满足perl脚本的要求，卡在最后生成这个可加载到tripal文件的这步

因此去把跟这个at vv文件相关的内容里面的都根据chr文件替换了一下，确保一致。再次重复上面的步骤，虽然仍然报错了，但是报的是那些命名不太规范的类似于17r这种染色体命名的问题，其他染色体命名都跑出结果了：

[root@f6f5e8d1a851 tripal_synview-master]# perl ./syntenyTool.pl -t mcscanx_tripal -a at -b vv -c 1 -d 2 /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.chr /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.gff /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.collinearity /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block > /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block.tripal.txt
[ERR]undef ref2 chr vv chr18r

五、前端页面测试模块

于是现在拿这个半成品去走tripal的流程

首先add content -> Synteny File

接着在**Human-readable Name for synteny File** *这部分里面写入一个类似 synteny analysis for at and vv 的内容

最后写入咱们上面跑出来的 at_vv.block.tripal.txt 的绝对文件路径

记得选genome，就是咱们前面对应的organism1和2

我这里还勾选了一下 Insert/update synteny files to chado database ，因为说明文档这里说了：Finally, do not forget check option for "Insert/update synteny files to chado database" and run tripal jobs to insert blocks and collinearity genes to datbase.试试：

跳出一个页面，需要drush一下，又错了

2024-05-11 14:32:39: Job ID 148.
2024-05-11 14:32:39: Calling: synview_parse_synfile(92, 92, /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block.tripal.txt)

NOTE: Checking the synteny file: 
  /var/www/html/sites/default/files/raw/MCScanX-master/data/at_vv.block.tripal.txt

FAILED: can not pass feature check 92 chr1 1 ...
Drush command terminated abnormally due to an unrecoverable error.

推测是物种id和这个at vv示例文件不一致导致的，到这基本感觉就跑通了，按照标准的要求的格式来整体做一遍应该就ok了，我后面再筛选两套比较标准的数据去从头到尾跑一边吧

建议适当的加一些图

如果有图了，就可以直观看到结果。