如何根据染色体坐标批量提取对应的DNA序列（bedtools）

这一篇小笔记是在我处理自己的数据的时候遇到的问题，经过查阅资料解决了，故记录下来。

比如现在：你需查找一段序列，比如说小鼠的chr10:105280000-105280550，我相信学生物的童鞋应该都知道应该怎么获得DNA序列，但是如果当我有上千条序列需要获得并把它们放在同一个fasta文件里的时候，应该怎么做呢？

方法如下：

Step1 你需要先拿到差异peaks

从ATAC-seq数据中分析得到的差异peaks，当然同样适用于chip-seq的数据：

> head(Q_T_Q_V_psig)
log2 fold change (MLE): condition Q_T vs Q_V 
Wald test p-value: condition Q_T vs Q_V 
DataFrame with 6 rows and 6 columns
                                  baseMean   log2FoldChange             lfcSE             stat               pvalue               padj
                                 <numeric>        <numeric>         <numeric>        <numeric>            <numeric>          <numeric>
chr10:105280000-105280550 9.37332336803852 3.53597545377286  1.04368265396652 3.38797951689085 0.000704095222812605                 NA
chr10:105469500-105469950 13.4771925765997 2.29247043028083  0.84387643836949  2.7165949018677  0.00659572827117144  0.353853732289383
chr10:107658700-107659600 58.6902716992164 1.13822876920294 0.433762615192843 2.62408222685789  0.00868828070349121  0.389236862222314
chr10:108153100-108153800 38.0887166659289 2.20816858694232 0.491153263586065 4.49588499284274 6.92811780142922e-06 0.0204186477574836
chr10:108452300-108452850 15.8518595117141 3.34668830781992 0.888088721665981 3.76841663020088 0.000164286335774994 0.0816701024146028
chr10:110183500-110184500 50.6887818370494 1.31842727975894 0.416836274255839 3.16293797153972  0.00156185605940921  0.204583310689789

随后把这个差异peaks表存成了csv文件。

Step2 在命令行里将一列分割成多列

从csv文件里提取某一列并存到另一个文件里，例如提取第一列（$1）:

$ awk -F"," '{print $1}' file.csv >> file1.csv

在linux里，再把csv文件的一列按照冒号分隔成两列，并存到另一个文件里：

$ sed 's/:/\t/' file.cvs > file1.csv

同样的，把csv文件的一列按照“-”短横杠进行分割：

$ sed 's/-/\t/' test1 > test2

处理后的文件：

$ head Q_T_Q_V_p01_chrlocation.csv
chr10   105280000       105280550
chr10   105469500       105469950
chr10   107658700       107659600
chr10   108153100       108153800
chr10   108452300       108452850
chr10   110183500       110184500
chr10   111636750       111637550
chr10   111684750       111685450
chr10   111840100       111840700
chr10   112000200       112000850

Step3 把csv改成bed后缀

直接鼠标右击“重命名”。
然后看一下bed文件：

$ head Q_T_Q_V_p01_chrlocation.bed 
chr10   105280000   105280550
chr10   105469500   105469950
chr10   107658700   107659600
chr10   108153100   108153800
chr10   108452300   108452850
chr10   110183500   110184500
chr10   111636750   111637550
chr10   111684750   111685450
chr10   111840100   111840700
chr10   112000200   112000850

Step4 安装bedtools

安装bedtools:

$ wget https://github.com/arq5x/bedtools2/releases/download/v2.29.1/bedtools-2.29.1.tar.gz
$ tar -zxvf bedtools-2.29.1.tar.gz
$ cd bedtools2
$ make

也有其他安装方法，见：Installation

Step5 根据染色体坐标位置批量提取相应序列

参考官网：getfasta
用bedtools根据已知的染色体坐标位置，获得fasta文件：

#-fi后面是你想从哪个fasta文件里提取序列，我这里是从小鼠基因组里提取
#-bed指的是你输入的文件类型是什么，这里我输入的是bed文件
#-bed后面跟的是上面我们拿到的染色体坐标文件
#-fo是输出结果储存的文件
$ bedtools getfasta -fi /media/yanfang/FYWD/RNA_seq/ref_genome/mm10.fa -bed Q_T_Q_V_p01_chrlocation.bed -fo Q_T_Q_V_p01_peakseq.fa

得到的fasta文件：