「 bedtools 」提取上游+gene+下游序列
1、bed文件格式介绍
BED文件每行至少包括chrom，chromStart，chromEnd三列必选;另外还可以添加额外的9列可选，这些列的顺序是固定的(之前一直以为时第五列，由于共线性里面分析的格式的第五列是正负，一直造成误解，啊啊啊啊啊)。
必选的三列：
1.chrom - 染色体的名称（例如chr3，chrY，chr2_random）或支架（例如scaffold10671）。
2.chromStart- 染色体或支架中特征的起始位置。染色体中的第一个碱基编号为0。
3.chromEnd- 染色体或支架中特征的结束位置。所述 chromEnd碱没有包括在特征的显示。例如，染色体的前100个碱基定义为chromStart = 0，chromEnd = 100，并跨越编号为0-99的碱基。
9个可选的BED字段：
4.name - 定义BED行的名称。当轨道打开到完全显示模式时，此标签显示在Genome浏览器窗口中BED行的左侧，或者在打包模式下直接显示在项目的左侧。
5.score - 得分在0到1000之间。如果此注释数据集的轨迹线useScore属性设置为1，则得分值将确定显示此要素的灰度级别（较高的数字=较深的灰色）。此表显示 Genome Browser将BED分数值转换为灰色阴影：
6.strand - 定义strand。要么“。” （=无绞线）或“+”或“ - ”。
7.thickStart- 绘制特征的起始位置（例如，基因显示中的起始密码子）。当没有厚部分时，thickStart和thickEnd通常设置为chromStart位置。
8.thickEnd - 绘制特征的结束位置（例如基因显示中的终止密码子）。
blockStart位置。此列表中的项目数应与blockCount相对应。
9.itemRgb- R，G，B形式的RGB值（例如255,0,0）。如果轨道行 itemRgb属性设置为“On”，则此RBG值将确定此BED行中包含的数据的显示颜色。注意：建议使用此属性的简单颜色方案（八种颜色或更少颜色），以避免压倒Genome浏览器和Internet浏览器的颜色资源。
10.blockCount- BED行中的块（外显子）数。
11.blockSizes- 块大小的逗号分隔列表。此列表中的项目数应与blockCount相对应。
12.blockStarts - 以逗号分隔的块开始列表。应该相对于chromStart计算所有blockStart**位置。此列表中的项目数应与blockCount相对应。
2.bedtools提取序列
bedtools：
[ Genome arithmetic ] #基因组区间运算
intersect-查看两个文件不同区间是否有 overlap；
closest-#找到和目标区间最近的特征区间；
coverage-计算特定区间的覆盖深度；
map-针对存在交叠区间的 B 文件的某一列应用函数如说求和、均值
genomecov-计算在整个基因组的覆盖深度；
merge-将重叠的或相邻的区间合并成一个区域；
cluster-类似 merge，但是只是增加一列标记，用于标明哪些行是聚集在一起的；
complement-提供一个区间，然后获得此区间与整个基因组不重叠的区域；
subtract-类似 complement，不是针对基因组，而是针对两个区域，去除掉一
个区域与另一个区域重叠部分；
slop-根据已有特征区间向外延展，可分别指定上下游延伸长度;
flank-根据现有区间，指定侧翼延伸长度，得到两侧翼位置的新区间，不包含现有区间；
sort-对区域进行排序；
random－根据 genome 文件随机生成指定长度指定个数的区域；
annotate-从其他多个文件中统计指定区间的覆盖深度；
2.1 提取gff文件的所有基因位置,并转换成bed格式
less ../01.data/EVM.final.gene.gff3|grep -w gene > genes.gff
convert2bed --input=gff --output=bed <genes.gff >genes.bed #调用的bedops，也可以自己用awk等提取
awk 示例 ：
less Mimosa_pudica_v1.gff |grep -w mRNA |awk '{print$1"\t"$4"\t"$5"\t"$9"\t"".""\t"$7}' |head
效果如下：
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 创建包含promoter位置的bed文件,up or down位置
slop-根据已有特征区间向外延展，可分别指定上下游延伸长度；
flank-根据现有区间，指定侧翼延伸长度，得到两侧翼位置的新区间，不包含现有区间；
一般默认启动子区域应该是上下游2kb，最大不超过5kb。
#up or down
#提取上游3k bed
bedtools flank -i genes.bed -g final.genome.fasta.len -l 3000 -r 0 -s > up_3k.promoters.bed
#提取gene bed
bedtools flank -i genes.bed -g final.genome.fasta.len -l 0 -r 2000 -s > down_2k.promoters.bed
#提取上下游的bed
bedtools flank -i genes.bed -g final.genome.fasta.len -l 3000 -r 2000 -s > up_down.promoters.bed
一起提取 up+gene+down序列，放在一条序列中，重点！！！我也不理解做启动子预测为什么要把基因序列加进去，但是有些人就喜欢这样，既然如此就满足各方需求。
bedtools slop -i genes.bed -g final.genome.fasta.len -l 3000 -r 2000 -s > up_3k.promoters.slop.bed
# -l 基因起始位置前多少bp
# # -r 基因后多少bp
# # -s
2.4 根据bed中的位置信息，在基因组序列中提取指定序列
分上下游提取;上下游一起提取，fa文件中id会一致;提取上游+gene+下游
#分上下游提取
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会一致
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序列
bedtools getfasta -s -fi ../01.data/03.Assembly_final/final.genome.fasta -bed genes.bed -fo J493.genes.fa -name
#提取上游+gene+下游，重点
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
如果超过坐标范围，bedtools会自己将坐标写到可提取到的坐标，如起始为0，末尾不足想要的长度时，提取的序列会用小写提示。
