這是唐海寶老師GitHub上的JCVI工具的非官方說(shuō)明書(shū)。
該工具集的功能非常多，但是教程資料目前看起來(lái)并不多，因此為了能讓更多人用上那么好用的工具，我就一邊探索，一邊寫(xiě)教程
這一篇文章教大家如何利用JCVI里面的工具繪制點(diǎn)圖，展現(xiàn)兩個(gè)物種之間的共線性關(guān)系。
在分析之前，你需要從PhytozomeV11 下載A.thaliana和Alyrata的CDS序列，保證文件夾里有如下內(nèi)容
Alyrata_384_v2.1.cds.fa.gz Athaliana_167_TAIR10.cds.fa.gz
Alyrata_384_v2.1.gene.gff3.gz Athaliana_167_TAIR10.gene.gff3.gz
準(zhǔn)備最長(zhǎng)CDS和BED文件
我們?cè)谧鯟DS相互比對(duì)的時(shí)候只需要有每個(gè)基因最長(zhǎng)的轉(zhuǎn)錄本即可，有兩種方法可以實(shí)現(xiàn)
方法1：自己寫(xiě)腳本
我用我寫(xiě)的一個(gè)腳本get_the_longest_transcripts.py提取每個(gè)基因的最長(zhǎng)轉(zhuǎn)錄本,見(jiàn) 基因組共線性工具M(jìn)CScanX使用說(shuō)明
zcat Alyrata_384_v2.1.gene.gff3.gz | python ~/scripts/python/get_the_longest_transcripts.py > aly_lst_gene.txt
zcat Athaliana_167_TAIR10.gene.gff3.gz | python ~/scripts/python/get_the_longest_transcripts.py > ath_lst_gene.txt
其中xxx_lst_gene.txt的格式如下, 第一列是基因名，第二列是mRNA編號(hào)，后面幾列是位置信息。
$ head ath_lst_gene.txt
AT4G19470.TAIR10 AT4G19470.1.TAIR10 Chr4 10612993 10614339 -
AT5G43860.TAIR10 AT5G43860.1.TAIR10 Chr5 17630450 17632312 +
AT1G68650.TAIR10 AT1G68650.1.TAIR10 Chr1 25775741 25777874 +
AT1G28050.TAIR10 AT1G28050.1.TAIR10 Chr1 9775528 9777810 -
AT3G59880.TAIR10 AT3G59880.1.TAIR10 Chr3 22120969 22121700 +
AT1G22030.TAIR10 AT1G22030.1.TAIR10 Chr1 7759164 7760556 -
AT5G24330.TAIR10 AT5G24330.1.TAIR10 Chr5 8295147 8297068 -
AT5G43990.TAIR10 AT5G43990.2.TAIR10 Chr5 17697889 17702005 +
AT1G11410.TAIR10 AT1G11410.1.TAIR10 Chr1 3841286 3844432 +
AT4G32890.TAIR10 AT4G32890.1.TAIR10 Chr4 15875470 15876762 +
由于基因名和mRNA編號(hào)里有在提取CDS不需要的內(nèi)容，因此要進(jìn)行刪除
sed -i 's/\.v2\.1//g' aly_lst_gene.txt
sed -i 's/\.TAIR10//g' ath_lst_gene.txt
之后我們就可以根據(jù)第二列進(jìn)行提取CDS
seqkit grep -f <(cut -f 2 ath_lst_gene.txt ) Athaliana_167_TAIR10.cds.fa.gz > ath.cds
seqkit grep -f <(cut -f 2 aly_lst_gene.txt ) Alyrata_384_v2.1.cds.fa.gz > aly.cds
提取的CDS編號(hào)里面也有一些不需要的內(nèi)容，所以也要?jiǎng)h除
sed -i 's/\.t.*//' aly.cds
sed -i 's/\..*//' ath.cds
此外還需要基因的位置信息的bed文件
awk '{print $3"\t"$4"\t"$5"\t"$1"\t0\t"$6}' ath_lst_gene.txt | sort -k4,4V > ath.bed
awk '{print $3"\t"$4"\t"$5"\t"$1"\t0\t"$6}' aly_lst_gene.txt | sort -k4,4V > aly.bed
基于JCVI工具集
當(dāng)然也可以參考「JCVI教程」如何基于編碼序列或蛋白序列進(jìn)行共線性分析來(lái)提取bed和cds序列，不需要用到我寫(xiě)的腳本。
python -m jcvi.formats.gff bed --type=mRNA --key=Name Athaliana_167_TAIR10.gene.gff3.gz > ath.bed
python -m jcvi.formats.gff bed --type=mRNA --key=Name Alyrata_384_v2.1.gene.gff3.gz > aly.bed
對(duì)bed文件中的基因進(jìn)行去重
python -m jcvi.formats.bed uniq ath.bed
python -m jcvi.formats.bed uniq aly.bed
這一步會(huì)得到aly.uniq.bed和ath.uniq.bed, 我們將其覆蓋原文件
mv ath.uniq.bed ath.bed
mv aly.uniq.bed aly.bed
根據(jù)bed文件提取cds里的序列
# Athaliana
seqkit grep -f <(cut -f 4 ath.bed ) Athaliana_167_TAIR10.cds.fa.gz | seqkit seq -i > ath.cds
# Alyrata
seqkit grep -f <(cut -f 4 aly.bed ) Alyrata_384_v2.1.cds.fa.gz | seqkit seq -i > aly.cds
無(wú)論是哪種方法，請(qǐng)保證最后有以下四個(gè)文件
$ ls ???.???
aly.bed aly.cds ath.bed ath.cds
BLAST比對(duì)
相對(duì)于上一步，這一步其實(shí)非常簡(jiǎn)單了
makeblastdb -in ath.cds -out db/ath -dbtype nucl
blastn -num_threads 20 -query aly.cds -db db/ath -outfmt 6 -evalue 1e-5 -num_alignments 5 > aly_ath.blast
用jcvi.compara.blastfilter對(duì)結(jié)果進(jìn)行過(guò)濾
python -m jcvi.compara.blastfilter --no_strip_names aly_ath.blast --sbed ath.bed --qbed aly.bed
運(yùn)行過(guò)程中有如下輸出信息
19:59:15 [base] Load file `aly.bed`
19:59:16 [base] Load file `ath.bed`
19:59:16 [blastfilter] Load BLAST file `aly_ath.blast` (total 49887 lines)
19:59:16 [base] Load file `aly_ath.blast`
19:59:16 [blastfilter] running the cscore filter (cscore>=0.70) ..
19:59:16 [blastfilter] after filter (42023->26531) ..
19:59:16 [blastfilter] running the local dups filter (tandem_Nmax=10) ..
19:59:16 [blastfilter] after filter (26531->24242) ..
最后輸出aly_ath.blast.filtered用于做圖
python -m jcvi.graphics.blastplot aly_ath.blast.filtered --sbed ath.bed --qbed aly.bed
最后點(diǎn)圖如下
