1、bed文件格式介紹
BED文件每行至少包括chrom，chromStart，chromEnd三列必選;另外還可以添加額外的9列可選浴井，這些列的順序是固定的(之前一直以為時第五列，由于共線性里面分析的格式的第五列是正負歧胁，一直造成誤解滋饲，啊啊啊啊啊)。
必選的三列：
1.chrom - 染色體的名稱（例如chr3喊巍，chrY，chr2_random）或支架（例如scaffold10671）箍鼓。
2.chromStart - 染色體或支架中特征的起始位置崭参。。
3.chromEnd - 染色體或支架中特征的結束位置款咖。所述 chromEnd堿沒有包括在特征的顯示何暮。例如，染色體的前100個堿基定義為chromStart = 0铐殃，chromEnd = 100海洼，并跨越編號為0-99的堿基。
9個可選的BED字段：
4.name - 定義BED行的名稱富腊。當軌道打開到完全顯示模式時坏逢，此標簽顯示在Genome瀏覽器窗口中BED行的左側，或者在打包模式下直接顯示在項目的左側。
5.score - 得分在0到1000之間是整。如果此注釋數(shù)據(jù)集的軌跡線useScore屬性設置為1肖揣，則得分值將確定顯示此要素的灰度級別（較高的數(shù)字=較深的灰色）。此表顯示 Genome Browser將BED分數(shù)值轉換為灰色陰影：
6.strand - 定義strand浮入。要么“塑悼⊙遥” （=無絞線）或“+”或“ - ”。
7.thickStart - 繪制特征的起始位置（例如，基因顯示中的起始密碼子）固以。當沒有厚部分時，thickStart和thickEnd通常設置為chromStart位置夭委。
8.thickEnd - 繪制特征的結束位置（例如基因顯示中的終止密碼子）堵泽。
blockStart位置。此列表中的項目數(shù)應與blockCount相對應赴蝇。
9.itemRgb - R菩浙，G，B形式的RGB值（例如255,0,0）句伶。如果軌道行 itemRgb屬性設置為“On”劲蜻，則此RBG值將確定此BED行中包含的數(shù)據(jù)的顯示顏色。注意：建議使用此屬性的簡單顏色方案（八種顏色或更少顏色）考余，以避免壓倒Genome瀏覽器和Internet瀏覽器的顏色資源先嬉。
10.blockCount - BED行中的塊（外顯子）數(shù)。
11.blockSizes - 塊大小的逗號分隔列表楚堤。此列表中的項目數(shù)應與blockCount相對應疫蔓。
12.blockStarts - 以逗號分隔的塊開始列表。應該相對于chromStart計算所有blockStart**位置身冬。此列表中的項目數(shù)應與blockCount相對應衅胀。
2.bedtools提取序列
bedtools：
[ Genome arithmetic ] #基因組區(qū)間運算
intersect-查看兩個文件不同區(qū)間是否有 overlap；
closest-#找到和目標區(qū)間最近的特征區(qū)間酥筝；
coverage-計算特定區(qū)間的覆蓋深度滚躯；
map-針對存在交疊區(qū)間的 B 文件的某一列應用函數(shù)如說求和、均值
genomecov-計算在整個基因組的覆蓋深度嘿歌；
merge-將重疊的或相鄰的區(qū)間合并成一個區(qū)域掸掏；
cluster-類似 merge，但是只是增加一列標記宙帝，用于標明哪些行是聚集在一起的丧凤；
complement-提供一個區(qū)間，然后獲得此區(qū)間與整個基因組不重疊的區(qū)域步脓；
subtract-類似 complement愿待，不是針對基因組浩螺，而是針對兩個區(qū)域，去除掉一
個區(qū)域與另一個區(qū)域重疊部分呼盆；
slop-根據(jù)已有特征區(qū)間向外延展年扩，可分別指定上下游延伸長度;
flank-根據(jù)現(xiàn)有區(qū)間，指定側翼延伸長度访圃，得到兩側翼位置的新區(qū)間厨幻，不包含現(xiàn)有區(qū)間；
sort-對區(qū)域進行排序腿时；
random－根據(jù) genome 文件隨機生成指定長度指定個數(shù)的區(qū)域况脆；
annotate-從其他多個文件中統(tǒng)計指定區(qū)間的覆蓋深度；
2.1 提取gff文件的所有基因位置,并轉換成bed格式
bedtools起始坐標是0批糟，注意$4-1格了，gff一般起始坐標是1。例如bedtolls寫0-1提取的是0位的堿基徽鼎，所以需要$4-1盛末，（一般我們寫0-1是想提取的是0和1位置的堿基）
convert2bed ，gff2bed 提取則直接是$4-1的結果否淤。
###將標準注釋gff3文件轉換得到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 #結果同上
###方法2
less EVM.final.gene.gff3 |grep -w gene|awk '{print$1"\t"$4-1"\t"$5"\t"$9"\t"".""\t"$7}' >genes.gff
：bedtools起始坐標是0石抡，注意$4-1檐嚣，gff一般起始坐標是1。例如bedtolls寫0-1提取的是0位的堿基啰扛，所以需要$4-1嚎京，（一般我們寫0-1是想提取的是0和1位置的堿基）
convert2bed 提取則直接是$4-1的結果。
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 計算染色體長度
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位置
slop-根據(jù)已有特征區(qū)間向外延展隐解，可分別指定上下游延伸長度鞍帝；
flank-根據(jù)現(xiàn)有區(qū)間，指定側翼延伸長度煞茫，得到兩側翼位置的新區(qū)間膜眠，不包含現(xiàn)有區(qū)間；
一般默認啟動子區(qū)域應該是上下游2kb溜嗜，最大不超過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序列炸宵，放在一條序列中，重點９瓤邸Ｍ寥捎琐！我也不理解做啟動子預測為什么要把基因序列加進去，但是有些人就喜歡這樣裹匙，既然如此就滿足各方需求瑞凑。
#提取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ù)說明
# -l 基因起始位置前多少bp
# -r 基因后多少bp
# -s 考慮正負鏈
2.4 根據(jù)bed中的位置信息，在基因組序列中提取指定序列
結果1：分上下游提取概页。
結果2：上下游一起提取籽御，fa文件中id會一致。
結果3：提取上游+gene+下游惰匙。
#分上下游提燃继汀：結果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會一致：結果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+下游项鬼，重點：結果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序列哑梳，不考慮延長基因坐標左右的邊位置。
bedtools getfasta -s -fi ../01.data/03.Assembly_final/final.genome.fasta -bed genes.bed -fo genes.fa -name
###參數(shù)說明：
-name 顯示名字绘盟，即bed的第四列的名字鸠真。不寫則顯示坐標的范圍
-s strand #考慮正負鏈條
-fi 基因組fa文件
-bed 準備好的bed格式文件
-fo 輸出文件名
2.5 包裝成shell
export PATH=/share/nas1/pengzw/software/bedops-v2.4.38/bin/:$PATH
#step0:數(shù)據(jù)準備
gff=Chr_genome_final_gene.gff3
genome=Lachesis_assembly_changed.fa
fai=Lachesis_assembly_changed.fa.fai
#step1:將標準注釋gff3文件轉換得到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 #結果同上
#step2準備基因組長度文件
#samtools faidx $genome >fai
cut -f 1,2 $fai >genome.len
len=genome.len
#step3:getfasta
#flank 分別提取上游撤嫩，下游序列，因為寫在儀器蠢终，id會一樣序攘，無法區(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
