這一篇小筆記是在我處理自己的數(shù)據(jù)的時候遇到的問題，經(jīng)過查閱資料解決了，故記錄下來。
比如現(xiàn)在：你需查找一段序列，比如說小鼠的chr10:105280000-105280550，我相信學(xué)生物的童鞋應(yīng)該都知道應(yīng)該怎么獲得DNA序列，但是如果當我有上千條序列需要獲得并把它們放在同一個fasta文件里的時候，應(yīng)該怎么做呢？
方法如下：
Step1 你需要先拿到差異peaks
從ATAC-seq數(shù)據(jù)中分析得到的差異peaks，當然同樣適用于chip-seq的數(shù)據(jù)：
> 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 根據(jù)染色體坐標位置批量提取相應(yīng)序列
參考官網(wǎng)：getfasta
用bedtools根據(jù)已知的染色體坐標位置，獲得fasta文件：
#-fi后面是你想從哪個fasta文件里提取序列，我這里是從小鼠基因組里提取
#-bed指的是你輸入的文件類型是什么，這里我輸入的是bed文件
#-bed后面跟的是上面我們拿到的染色體坐標文件
#-fo是輸出結(jié)果儲存的文件
$ 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文件：
