用BEDtools/Python序列截取

在用MEME-chip进行模体搜索时，需要将fasta格式的文件上传，而上传的数据又必须是大于50bp且等长的序列，这个时候就要对数据进行序列截取。此处将介绍两种截取序列的方法。

一、BEDtools截取序列

Bedtools是一款可以对genomic features进行比较、相关操作和注释的工具，目前版本已经有三十多个工具／命令用以实现各种不同的功能，可以针对bed、vcf、gff等格式的文件进行处理。

BEDtools的官方主页地址为https://bedtools.readthedocs.io/en/latest/content/bedtools-suite.html
在bedtools中众多功能中，getfasta可截取序列并获得fasta文件。

getfasta
进行序列截取时应该已经有了全基因组序列数据和peak数据，首先用excle表格对数据进行初步处理，以Aoaf细胞系结合CTCF的peak数据为例
数据示例
若想要截取包含peak，100bp长的序列，首先用公式取start和end的中点=ROUND(AVERAGE(start1,end1),0)
取中点
然后再分别计算中点减50和中点加五十，这样我们就得到了100bp区间的包含peak的数据
100bp区间 然后我们将A，L，M列保存为一个新的bed文件。接下来我们打开bedtools开始序列截取。
首先将24条染色体的序列信息下载好，然后使用命令

bedtools getfasta -fi /usr/chailu/2020job/motif_discovery/jiyinzu/chr1.fa  -bed /usr/chailu/2020job/motif_discovery/shuru/Gm12892.bed  -fo /usr/chailu/2020job/motif_discovery/shuchu/1.fa

以上命令只截取了chr1的，全部截取需要更改文件名字进行下一次截取。
bedtools getfasta意思是使用bedtools下getfasta功能。
-fi意思是输入的染色体序列信息的文件，后面跟着的是文件的位置信息。
-bed意思是输入的bed文件，跟着的是文件位置信息。
-fo意思是输出文件所在位置及对输出文件命名。

最终我们会得到24个fasta文件分别对应各个染色体的截取结果，最终我们用一个简单的python程序将这些文件合并，得到总的fasta文件，将这个文件上传到MEME-chip上就可以进行接下来的模体搜索啦。

二、Python程序截取序列

考虑到有时候Linux系统的复杂性，或者有时候服务器不好使了，我们的备选就是人生苦短的Python啦，前提也是首先得有全基因组序列信息，和peak文件，对peak文件进行处理同上，命名为end.txt

end.txt

然后用以下python程序进行序列截取

fp=open(r'F:\2020myjob\jiyinzu\shuru\end.txt','rt') #位置信息文件
ff=open(r'F:\2020myjob\jiyinzu\chr1.fa','rt')  #多行变一行后的序列文件导入
fw=open(r'F:\2020myjob\jiyinzu\shuchu/1.txt','wt')   #改对应染色体名字，输出文件，result
line_fasta=ff.readlines()
line_position=fp.readlines()
seq=[]
sequence=''
for line_a in line_fasta:
    line_a=line_a.strip().replace('|',' ')
    seq.append(line_a)
sequence=sequence.join(seq)

result=[]
for line in line_position:
    line = line.strip().split('\t')
    if 'chr1' in line:     #染色体条数
        start=line[1]#待截取序列起始位点
        end=line[2]#待截取序列终止位点
        fragment=sequence[int(start):int(end)]
        print(fragment,'\t',file=fw)

fp.close()
ff.close()
fw.close()

print('yeah')

最终我们会得到截取出来的序列seq信息，但是并不是fasta格式的文件，首先我们将这24个序列文件进行合并，合并文件的python程序为

import os
#获取目标文件夹的路径
os.getcwd()
os.chdir('F:\\2020myjob\\jiyinzu')
meragefiledir = os.getcwd()+'\\shuchu'
#获取当前文件夹中的文件名称列表
filenames=os.listdir(meragefiledir)
file=open('Hmf.txt','w')
#先遍历文件名
for filename in filenames:
    filepath=meragefiledir+'\\'+filename
    #遍历单个文件，读取行数
    for line in open(filepath):
        file.writelines(line)
    file.write('\n') 
file.close()

这个程序运行完之后会产生空行，接下来再用去除空行程序来将中间的空行去除

def delblankline(infile, outfile):
 infopen = open(infile, 'r',encoding="utf-8")
 outfopen = open(outfile, 'w',encoding="utf-8")
 
 lines = infopen.readlines()
 for line in lines:
  if line.split():
   outfopen.writelines(line)
  else:
   outfopen.writelines("")
 
 infopen.close()
 outfopen.close()
 
delblankline("Hmf.txt", "Hmf_seq.txt")#Hmf.txt是没去空行前的文件，Hmf_seq.txt是去除空行后的文件

此处我们已经得到了截取出来的序列信息，接下来我们再用excle对之前的end.txt再进行一次处理输入公式=A1&":"&B1&"-"&C1来得到编号文件。

A.txt
将D列再保存为一个txt文件，命名为A.txt，此时我们就可用seq序列文件和编号文件A.txt,生成一个fasta格式的文件，同样用python程序来达到这个目的

import os
os.getcwd()
os.chdir('F:\\')
f_ = open('F:/Hee_seq.txt','r')#序列文件
n=0
list1=[]
for i in f_.readlines():
    n+=1
    s=i.strip()
    list1.append(s)
f_.close()
 
ff_ = open('F:/A.txt','r')#编号文件
m=0
list2=[]
for i in ff_.readlines():
    m+=1
    s=i.strip()
    list2.append(s)
ff_.close()
 
fff_=open('F:/Hee.fa','w')
for i in range(n):
    s='>' + list2[i]+'\n'+list1[i]
    fff_.write(s+'\n')
    #print(s)
fff_.close()    
print('OK!')

得到的fasta文件就可以进行模体搜索了。

Tips：切记对初始的peak排序把染色体排序成按顺序1-22，X，Y，因为
Python内置的排序方式为1，10，11，12，13，14，15，16，17，18，19，2，20，21，22，3，4，5，6，7，8，9，X，Y，所以合并文件的时候如果你的命名是数字的话就肯定会是这样的排序，合成fasta格式文件的时候就会产生编号和序列并不匹配的错误，若要避免这个情况的发生，可以在对输出文件命名时以字母开头，按abcd···的顺序输出，或者不对peak文件排序就可避免这个错误啦。

注：文中所有python程序均来自网络或师姐的编辑