一篇不知道拖了多久的筆記……欠下的債早晚還是要還的……
0.軟件下載及安裝
官網(wǎng):
參考教程:
其實(shí)之前已經(jīng)有兩位前輩寫過教程了的秀撇,我這只是再給需要的朋友舉個(gè)例子打個(gè)樣。詳情請(qǐng)戳下面兩個(gè)鏈接
下載軟件:
解壓:
unzip MCScanX.zip
cd到MCScanX后捌袜,按照
介紹的方法虏等,給
msa.h
dissect_multiple_alignment.h
detect_collinear_tandem_arrays.h
這三個(gè)文件前面添加上#include <unistd.h>
即弄唧,用vim將文件打開后在開頭加上上述的參數(shù)段。這個(gè)錯(cuò)誤的原因是霍衫,MCScanX 不支持64位系統(tǒng)候引。
然后make一下就安裝好了
1.文件準(zhǔn)備
從TAIR上下載TAIR10的pep和gff3數(shù)據(jù),要先用notepad++的正則替換功能(vim的批量替換如果會(huì)用的話也是可以的)將蛋白質(zhì)的基因ID后面的注釋信息給刪掉敦跌,否則會(huì)影響后續(xù)的工作澄干。接下來處理gff3文件。
先看一眼原始的數(shù)據(jù)長啥樣:
我們需要的是第1柠傍、4麸俘、5列和第9列的ID部分
用awk、grep和vim的全局替換功能可以得到所需要的簡化版gff3文件:
cat TAIR10.gff3 |awk '/gene/ {print $1 "\t" $3 "\t" $4 "\t" $5 "\t" $9 }'> tmp.gff3
$3是用來標(biāo)記用的惧笛,我們只取含有g(shù)ene的列
但是會(huì)有一些問題从媚,比如:
在搜索gene的時(shí)候只要包含gene的詞就會(huì)被檢索到,于是有以下三個(gè)混進(jìn)來的:
pseudogene
transposable_element_gene
mRNA_TE_gene
于是我就:
grep -v 'mRNA_TE_gene' tmp2.gff3 > tmp3.gff3
grep -v 'pseudogene' tmp3.gff3 > tmp4.gff3
grep -v 'transposable_element_gene' tmp4.gff3 > tmp5.gff3
由于當(dāng)時(shí)學(xué)藝不精比較菜就只能先暫時(shí)這樣用一下了……
現(xiàn)在想想好像在上一步awk的時(shí)候可以:
cat TAIR10.gff3 |awk '/^gene/ {print $2 "\t" $0}' > tmp.gff3
加上一個(gè)“^”表示以gene開頭拜效,應(yīng)該就可以免去上面那繁瑣的幾步了,最后生成的時(shí)候再倒一下回來就好了。
可以用
grep -v "^gene" tmp3.gff3
來檢測是否還包含有非gene的行各谚。
經(jīng)過一系列的awk操作之后獲得最終的結(jié)果:
另一個(gè)物種如法炮制即可紧憾。
updates:也可以這么寫,三步搞定:
awk '$3 == "gene" {print $1 "\t" $3 "\t" $4 "\t" $5 "\t" $9}' TAIR10.gff3 > new.gff3
用vim全局替換處理ID列的冗余信息
用awk將信息位置調(diào)換正確
2.建庫
makeblastdb -in AT.pep -dbtype prot -parse_seqids -out ATdb
因?yàn)槲倚枰龅氖莔aca對(duì)擬南芥的共線性比較昌渤,所以只需要對(duì)擬南芥建庫即可赴穗。(如果需要物種內(nèi)也進(jìn)行共線性比較的話這里需要把兩個(gè)物種的蛋白質(zhì)文件cat到一起來建庫)
3.比對(duì)
nohup blastp -query maca.pep -db ATdb -out AT_maca.blast -evalue 1e-10 -num_threads 30 -outfmt 6 -num_alignments 5 &
參數(shù)稍微解釋一下。
-query 被比對(duì)的物種的蛋白質(zhì)文件
-db 上一步建的庫
-out 輸出結(jié)果
-evalue 設(shè)置輸出結(jié)果的e-value值愈涩,值越小相似度越高望抽。
-num_threads 比對(duì)所用的線程數(shù)
-outfmt 輸出的格式
-num_alignments 輸出比對(duì)上的序列的最大值條目數(shù)
這一步還是要點(diǎn)時(shí)間的。記得掛nohup履婉。
4.運(yùn)行MCScanX
將之前的兩個(gè)物種的gff3文件cat成一個(gè)文件后命名為AT_maca.gff
將上一步生成的AT_maca.blast和AT_maca.gff放到新建的AT_maca文件夾內(nèi)
MCScanX ./AT_maca
之后就會(huì)生成一個(gè)叫AT_maca.collinearity的文件和一個(gè)AT_maca.html文件夾煤篙。
AT_maca.html里面是這樣滴:
5.下游分析
java circle_plotter -g AT_maca.gff -s AT_maca.collinearity -c circle.ctl -o AT_maca_circle.png
這一步需要修改修改circle.ctl文件
只要把你希望畫出來的序列的ID給寫進(jìn)去就好了
可以畫dot圖,bar圖和circle圖毁腿,個(gè)人覺得circle圖最好看了~
大概就是這樣了辑奈。如果有什么錯(cuò)誤的地方歡迎討論指出喲。