「 bedtools 」提取上游+gene+下游序列

1、bed文件格式介紹

BED文件每行至少包括chrom，chromStart，chromEnd三列必選;另外還可以添加額外的9列可選，這些列的順序是固定的(之前一直以為時(shí)第五列，由于共線性里面分析的格式的第五列是正負(fù)，一直造成誤解，啊啊啊啊啊)。

必選的三列：
1.chrom - 染色體的名稱（例如chr3，chrY，chr2_random）或支架（例如scaffold10671）。
2.chromStart - 染色體或支架中特征的起始位置。\color{red}{染色體中的第一個(gè)堿基編號(hào)為0}。
3.chromEnd - 染色體或支架中特征的結(jié)束位置。所述 chromEnd堿沒(méi)有包括在特征的顯示。例如，染色體的前100個(gè)堿基定義為chromStart = 0，chromEnd = 100，并跨越編號(hào)為0-99的堿基。
9個(gè)可選的BED字段：
4.name - 定義BED行的名稱。當(dāng)軌道打開到完全顯示模式時(shí)，此標(biāo)簽顯示在Genome瀏覽器窗口中BED行的左側(cè)，或者在打包模式下直接顯示在項(xiàng)目的左側(cè)。
5.score - 得分在0到1000之間。如果此注釋數(shù)據(jù)集的軌跡線useScore屬性設(shè)置為1，則得分值將確定顯示此要素的灰度級(jí)別（較高的數(shù)字=較深的灰色）。此表顯示 Genome Browser將BED分?jǐn)?shù)值轉(zhuǎn)換為灰色陰影：
6.strand - 定義strand。要么“?！?（=無(wú)絞線）或“+”或“ - ”。
7.thickStart - 繪制特征的起始位置（例如，基因顯示中的起始密碼子）。當(dāng)沒(méi)有厚部分時(shí)，thickStart和thickEnd通常設(shè)置為chromStart位置。
8.thickEnd - 繪制特征的結(jié)束位置（例如基因顯示中的終止密碼子）。
blockStart位置。此列表中的項(xiàng)目數(shù)應(yīng)與blockCount相對(duì)應(yīng)。
9.itemRgb - R，G，B形式的RGB值（例如255,0,0）。如果軌道行 itemRgb屬性設(shè)置為“On”，則此RBG值將確定此BED行中包含的數(shù)據(jù)的顯示顏色。注意：建議使用此屬性的簡(jiǎn)單顏色方案（八種顏色或更少顏色），以避免壓倒Genome瀏覽器和Internet瀏覽器的顏色資源。
10.blockCount - BED行中的塊（外顯子）數(shù)。
11.blockSizes - 塊大小的逗號(hào)分隔列表。此列表中的項(xiàng)目數(shù)應(yīng)與blockCount相對(duì)應(yīng)。
12.blockStarts - 以逗號(hào)分隔的塊開始列表。應(yīng)該相對(duì)于chromStart計(jì)算所有blockStart**位置。此列表中的項(xiàng)目數(shù)應(yīng)與blockCount相對(duì)應(yīng)。

2.bedtools提取序列

bedtools：
[ Genome arithmetic ] #基因組區(qū)間運(yùn)算
intersect-查看兩個(gè)文件不同區(qū)間是否有 overlap；
closest-#找到和目標(biāo)區(qū)間最近的特征區(qū)間；
coverage-計(jì)算特定區(qū)間的覆蓋深度；
map-針對(duì)存在交疊區(qū)間的 B 文件的某一列應(yīng)用函數(shù)如說(shuō)求和、均值
genomecov-計(jì)算在整個(gè)基因組的覆蓋深度；
merge-將重疊的或相鄰的區(qū)間合并成一個(gè)區(qū)域；
cluster-類似 merge，但是只是增加一列標(biāo)記，用于標(biāo)明哪些行是聚集在一起的；
complement-提供一個(gè)區(qū)間，然后獲得此區(qū)間與整個(gè)基因組不重疊的區(qū)域；
subtract-類似 complement，不是針對(duì)基因組，而是針對(duì)兩個(gè)區(qū)域，去除掉一
個(gè)區(qū)域與另一個(gè)區(qū)域重疊部分；
slop-根據(jù)已有特征區(qū)間向外延展，可分別指定上下游延伸長(zhǎng)度;
flank-根據(jù)現(xiàn)有區(qū)間，指定側(cè)翼延伸長(zhǎng)度，得到兩側(cè)翼位置的新區(qū)間，不包含現(xiàn)有區(qū)間；
sort-對(duì)區(qū)域進(jìn)行排序；
random－根據(jù) genome 文件隨機(jī)生成指定長(zhǎng)度指定個(gè)數(shù)的區(qū)域；
annotate-從其他多個(gè)文件中統(tǒng)計(jì)指定區(qū)間的覆蓋深度；

2.1 提取gff文件的所有基因位置,并轉(zhuǎn)換成bed格式

\color{red}{注意事項(xiàng)}
bedtools起始坐標(biāo)是0，注意$4-1，gff一般起始坐標(biāo)是1。例如bedtolls寫0-1提取的是0位的堿基，所以需要$4-1，（一般我們寫0-1是想提取的是0和1位置的堿基）
convert2bed ，gff2bed 提取則直接是$4-1的結(jié)果。


###將標(biāo)準(zhǔn)注釋gff3文件轉(zhuǎn)換得到bed格式的gff文件
###方法1  #調(diào)用的bedops，也可以自己
awk '{if($3~/^gene$/)print}' EVM.gff3 > genes.gff && convert2bed  --input=gff --output=bed  < gene.gff >genes.bed #調(diào)用的bedops，也可以自己用awk
awk '{if($3~/^gene$/)print}' EVM.gff3 > genes.gff && gff2bed <genes.gff> genes.bed #結(jié)果同上
###方法2
less EVM.final.gene.gff3 |grep -w gene|awk '{print$1"\t"$4-1"\t"$5"\t"$9"\t"".""\t"$7}'  >genes.gff

\color{red}{注意事項(xiàng)}：bedtools起始坐標(biāo)是0，注意$4-1，gff一般起始坐標(biāo)是1。例如bedtolls寫0-1提取的是0位的堿基，所以需要$4-1，（一般我們寫0-1是想提取的是0和1位置的堿基）
convert2bed 提取則直接是$4-1的結(jié)果。

bed格式效果如下：
scaffold10018_cov214    9657    10808   ID=Mimpu10018S00001     .  -
scaffold1001_cov228     131443  132576  ID=Mimpu1001S00002      .  +
scaffold1001_cov228     134493  139136  ID=Mimpu1001S00003      .  -
scaffold1001_cov228     38129   38701   ID=Mimpu1001S00004      .  +
scaffold10020_cov190    52371   53006   ID=Mimpu10020S00005     .  +
scaffold10020_cov190    72945   74067   ID=Mimpu10020S00006     .  -
scaffold10020_cov190    171722  173389  ID=Mimpu10020S00007     .  +
scaffold10020_cov190    137289  137993  ID=Mimpu10020S00008     .  -
2.2 計(jì)算染色體長(zhǎng)度
samtools faidx ../01.data/03.Assembly_final/final.genome.fasta
cut -f 1,2 ../01.data/03.Assembly_final/final.genome.fasta.fai >final.genome.fasta.len
2.3 創(chuàng)建包含promoter位置的bed文件,up or down位置

\color{red}{注意事項(xiàng)}
slop-根據(jù)已有特征區(qū)間向外延展，可分別指定上下游延伸長(zhǎng)度；
flank-根據(jù)現(xiàn)有區(qū)間，指定側(cè)翼延伸長(zhǎng)度，得到兩側(cè)翼位置的新區(qū)間，不包含現(xiàn)有區(qū)間；
一般默認(rèn)啟動(dòng)子區(qū)域應(yīng)該是上下游2kb，最大不超過(guò)5kb。

#up or down 
# 提取基因上游3k 區(qū)間
bedtools flank -i genes.bed -g final.genome.fasta.len  -l 3000 -r 0 -s > up_3k.promoters.bed
# 提取基因下游2k的區(qū)間
bedtools flank -i genes.bed -g final.genome.fasta.len  -l 0 -r 2000 -s > down_2k.promoters.bed
# 提取基因上游3k，下游2k區(qū)間
bedtools flank -i genes.bed -g final.genome.fasta.len  -l 3000 -r 2000 -s > up_down.promoters.bed

一起提取 up+gene+down序列，放在一條序列中，重點(diǎn)！??！我也不理解做啟動(dòng)子預(yù)測(cè)為什么要把基因序列加進(jìn)去，但是有些人就喜歡這樣，既然如此就滿足各方需求。

#提取up+gene+down區(qū)間
bedtools slop  -i genes.bed -g final.genome.fasta.len  -l 3000 -r 2000 -s > up_3k.promoters.slop.bed
#### 參數(shù)說(shuō)明
# -l 基因起始位置前多少bp
# -r 基因后多少bp
# -s 考慮正負(fù)鏈
2.4 根據(jù)bed中的位置信息，在基因組序列中提取指定序列

結(jié)果1：分上下游提取。
結(jié)果2：上下游一起提取，fa文件中id會(huì)一致。
結(jié)果3：提取上游+gene+下游。

#分上下游提?。航Y(jié)果1示例cmd
bedtools getfasta -s -fi ../01.data/03.Assembly_final/final.genome.fasta -bed up_3k.promoters.bed -fo up_3k.promoters.fa -name
bedtools getfasta -s -fi ../01.data/03.Assembly_final/final.genome.fasta -bed down_2k.promoters.bed -fo down_2k.promoters.fa -name
#上下游一起提取，fa文件中id會(huì)一致：結(jié)果2示例cmd
bedtools getfasta -s -fi ../01.data/03.Assembly_final/final.genome.fasta -bed up_down_2k.promoters.bed -fo down_2k.promoters.fa -name
#提取上游+gene+下游，重點(diǎn)：結(jié)果3示例cmd
bedtools getfasta -s -fi ../01.data/03.Assembly_final/final.genome.fasta -bed up_3k.promoters.slop.bed -fo J493.up_genes_down.fa -name

#提取gene序列，不考慮延長(zhǎng)基因坐標(biāo)左右的邊位置。
bedtools getfasta -s -fi ../01.data/03.Assembly_final/final.genome.fasta -bed genes.bed -fo genes.fa -name

###參數(shù)說(shuō)明：
-name 顯示名字，即bed的第四列的名字。不寫則顯示坐標(biāo)的范圍
-s strand #考慮正負(fù)鏈條
-fi 基因組fa文件
-bed 準(zhǔn)備好的bed格式文件
-fo 輸出文件名

\color{red}{注意事項(xiàng)}：如果超過(guò)坐標(biāo)范圍，bedtools會(huì)自己將坐標(biāo)寫到可提取到的坐標(biāo)，如起始為0，末尾不足想要的長(zhǎng)度時(shí)，提取的序列會(huì)用小寫提示。

2.5 包裝成shell
export PATH=/share/nas1/pengzw/software/bedops-v2.4.38/bin/:$PATH
#step0:數(shù)據(jù)準(zhǔn)備
gff=Chr_genome_final_gene.gff3
genome=Lachesis_assembly_changed.fa
fai=Lachesis_assembly_changed.fa.fai
#step1:將標(biāo)準(zhǔn)注釋gff3文件轉(zhuǎn)換得到bed格式的gff文件
convert2bed  --input=gff --output=bed  < $gff >all.bed && awk '{if($8~/^gene$/)print}' all.bed  >genes.bed
#awk '{if($3~/^gene$/)print}' $gff > genes.gff && gff2bed <genes.gff> genes.bed #結(jié)果同上

#step2準(zhǔn)備基因組長(zhǎng)度文件
#samtools faidx $genome >fai
cut -f 1,2 $fai >genome.len
len=genome.len

#step3:getfasta
#flank 分別提取上游，下游序列，因?yàn)閷懺趦x器，id會(huì)一樣，無(wú)法區(qū)分上下游
l=2000
r=2000
bedtools flank  -i genes.bed -g $len -l $l -r 0 -s > up.flank.bed
bedtools getfasta -s -fi $genome -bed up.flank.bed -fo up.gene.flank.promoter.fa -name 

bedtools flank  -i genes.bed -g $len -l 0 -r $r -s > down.flank.bed
bedtools getfasta -s -fi $genome -bed down.flank.bed -fo down.gene.flank.promoter.fa -name 
#slop 上下游延伸提取
l=2000
r=2000
bedtools slop  -i genes.bed -g $len -l $l -r $r -s > slop.bed
bedtools getfasta -s -fi $genome -bed slop.bed -fo all.gene.slop.promoter.fa -name

參考：
http://www.itdecent.cn/p/9208c3b89e44

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容