小白我最近做的事有些雜亂，是時(shí)候整理一波了诗祸。

序列的比對(duì)及提取已知ID的序列

故事是這樣開始的：

故事的主要內(nèi)容

如圖所示镀虐，我們遇到了三個(gè)問題屹培，當(dāng)然默穴，這么difficult的問題我是搞不定的啦，在謝大佬的幫助下這三個(gè)問題是已經(jīng)解決了褪秀，但我們不滿足于此蓄诽，我們想要把和這些探針序列匹配上的轉(zhuǎn)座子序列提取出來，以便我們進(jìn)一步設(shè)計(jì)更多的探針媒吗。謝大佬說剩下的工作比較簡單若专，就把他前面的筆記給我們，讓我們自己把轉(zhuǎn)座子序列提取出來蝴猪。于是乎调衰，照貓畫虎模式開啟。

1st Step

去NCBI把blast+下載到服務(wù)器賬號(hào)自阱，ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/嚎莉，在這個(gè)鏈接里選擇合適的版本，然后

wget -c 網(wǎng)址???????? —— 解壓 —— 添加環(huán)境變量 —— 搞定

注意：-c 斷點(diǎn)續(xù)傳 顧名思義沛豌，劃重點(diǎn)哦

此外還要在citrus.hzau.edu.cn下載甜橙csi.chromosome.fa序列, 并準(zhǔn)備探針序列的fasta格式文件趋箩。

2nd Step

將探針序列blast到甜橙基因組上赃额，生成格式7，以名為“CL”的探針為例叫确，格式7包含以下信息

格式7包含的信息

blast過程也很簡單跳芳，可以寫個(gè)腳本，提交到隊(duì)列中運(yùn)行

blast 腳本

第一行竹勉，打點(diǎn)當(dāng)前目錄下文件飞盆，并將其作為運(yùn)行參數(shù)

第二行，建立基因組的索引文件次乓，數(shù)據(jù)為Nt即DNA

這個(gè) “解析序列ID” 我還沒懂是神馬意思

第三行吓歇，將query序列blast到基因組上，輸出.out7文件

num_threads參數(shù)表明運(yùn)行blast所用CPU數(shù)票腰，應(yīng)該和向服務(wù)器申請(qǐng)的CPU數(shù)一致

3rd Step

從fish_csi.out7文件中提取出各個(gè)探針I(yè)dentity>80%, 且coverage>80%的blast結(jié)果

提取blast結(jié)果城看，寫入*_Ident0.8Cover0.8文件

這是以CL探針為例，我們就是這么一條探針一條探針地提取的杏慰，寫到腳本文件中運(yùn)行测柠，也可能有其他更快捷的辦法，但因?yàn)樘结様?shù)量不多缘滥，為了節(jié)約時(shí)間成本就沒再深入追究鹃愤。

知識(shí)點(diǎn)：對(duì)于“&&”的解釋可以參看下面的鏈接，意為 &&左邊的command1執(zhí)行成功(返回0表示成功)后完域，&&右邊的command2才能被執(zhí)行。與之對(duì)應(yīng)的還有“||”和“()”瘩将。

https://www.cnblogs.com/chenggang816/p/10303508.html

4th Step

將上一步生成的每條序列的篩選后的結(jié)果文件轉(zhuǎn)換成bed格式吟税，用于下步分析, 通過vim TransToBed.sh新建腳本文件，并在TransToBed.sh腳本文件中寫入

TransToBed.sh腳本文件中寫入的內(nèi)容

知識(shí)點(diǎn)：①Shell腳本for循環(huán)結(jié)構(gòu)如圖所示姿现；②“$”為申請(qǐng)參數(shù)肠仪，形成隊(duì)列；③“echo”為在Shell中打印备典，運(yùn)行腳本后可以看到終端打印出一行行*.Ident0.8Cover0.8的文件异旧，起到監(jiān)視識(shí)別文件是否正確的作用，此外“echo”輸出的結(jié)尾自帶換行符提佣，所以該命令結(jié)束后的 [賬戶名 目錄名]$ 是新開一行的吮蛹，而如果用“cat”命令顯示一個(gè)結(jié)尾無換行符的文本文件后 [賬戶名 目錄名]$ 是緊跟在文檔最后一個(gè)字符后面的，而不是新開一行拌屏，這在有利于在合并FASTA文件之判斷合并前的FASTA文件末尾是否有換行符潮针；④bed文件的分隔符為“\t”；⑤awk工具的If,else語句如圖倚喂；⑥圖中提供了將目錄下后綴相同的文件全部執(zhí)行操作后分別輸出到加了新后綴名的文件中每篷；⑦在目錄下對(duì)文件進(jìn)行批量操作時(shí)同一批操作的文件使用相同的后綴名，方便進(jìn)行批處理（我算是明白為何在NCBI的Gbrowse上檢索下載的序列文件有那么長而且整齊的文件名了）

5th Step

將轉(zhuǎn)座子的.gff3文件轉(zhuǎn)為bed格式。謝大佬不知從哪里拿來了citrus中的轉(zhuǎn)座子的注釋文件焦读，先看看長什么樣子吧

轉(zhuǎn)座子.gff3文件格式

由此看在citrus基因組不同位置上分布的轉(zhuǎn)座子是有不同的ID號(hào)的子库。

知道長啥樣，轉(zhuǎn)起格式來豈不是小菜一碟矗晃？

.gff3轉(zhuǎn)為.bed文件

6th Step

利用bedtools將每條序列的bed格式文件跟csi的TE的gff注釋文件的bed格式進(jìn)行intersect仑嗅，并判斷每條序列完全被某個(gè)TE包含(包含度為80%，填入的數(shù)值為已經(jīng)通過探針長度*0.8計(jì)算過的)

bedTools intersect是求兩個(gè)文件所表示的染色體區(qū)域的交集

bedtools 更多使用方式可參看以下大佬的鏈接

http://www.reibang.com/p/6c3b87301491

提取求交集后得到的轉(zhuǎn)座子的ID

7th Step

我們有了轉(zhuǎn)座子的ID名字下一步就是提取該轉(zhuǎn)座子在染色體上的位置信息了（雖然后來我發(fā)現(xiàn)其實(shí)用bedtools intersect 的 -wo選項(xiàng)輸出的內(nèi)容中有轉(zhuǎn)座子的位置信息喧兄，唉无畔，玩游戲還是要看攻略呀！）康康我是咋弄的吧吠冤。

剛開始浑彰，感覺自己還很機(jī)智

剛開始的錯(cuò)誤腳本

結(jié)果啥也沒輸出來，sort排過續(xù)后的文件根本沒有重復(fù)的拯辙，后來在uniq后面加了 -f 選項(xiàng)郭变，忽略幾行也不奏效，感覺可能是文件太大涯保，于是乎將文件拆分了诉濒，把各個(gè)探針對(duì)應(yīng)的轉(zhuǎn)座子ID分別提出來

把探針CL對(duì)應(yīng)的轉(zhuǎn)座子ID提出來

這里在第一列的名稱行加了“z”，是因?yàn)橛械奶结樏值腁SCII碼是小于chr的夕春，在sort的時(shí)候會(huì)被排在chr前面未荒，提取重復(fù)的時(shí)候就不能提取到位置信息的行

再把所有能和探針有交集的轉(zhuǎn)座子的注釋信息提出來，我借鑒了該鏈接下的回答及志，得到all_te.bed文件

https://segmentfault.com/q/1010000010766902

這個(gè)鏈接讓我學(xué)到了新的東西：如何用awk工具對(duì)兩個(gè)文件進(jìn)行不同卻有關(guān)聯(lián)的操作片排，感謝大佬。結(jié)果出來后發(fā)現(xiàn)得到的行數(shù)比repeat_te.bed的行數(shù)少一百多行速侈，我認(rèn)為應(yīng)該是有的轉(zhuǎn)座子中包含了兩個(gè)探針率寡，在awk雙文件操作時(shí)重復(fù)的下標(biāo)識(shí)別的是同一個(gè)參數(shù)，相當(dāng)于把ID一樣的轉(zhuǎn)座子又合并回同一個(gè)轉(zhuǎn)座子了倚搬。

之后就是重復(fù)第一步了

得到CL探針?biāo)鶎?duì)應(yīng)的位置信息

得到這些我就可以去基因組中提序列了

8th Step

首先冶共，將.loca文件轉(zhuǎn)為了.bed文件，之后用bedtools getfasta工具提取序列

提取序列

然后就可以用BioEdit看啦每界，雖然結(jié)果是 并沒有什么卵用 捅僵！ 但還是學(xué)到了點(diǎn)東西的。