去除测序reads中的接头：adaptor

之前用c写过一个程序，查找reads中是否包含了adaptor，如果检测到的话就过滤掉含有adaptor的reads，这次在过滤完数据之后发现接头序列比较多，为了提升组装效果，又不能很大地影响数据量，需要对接头进行截断处理，并过滤过短的reads，用python写了一个简短的程序，指定超过3个错配以内的匹配都认为匹配到，并且长度小于50bp的reads过滤，在以下程序基础上添加传入参数，可以适用比较多的情况（单端、双端、含有single等）：

import sys
import re
from Bio import SeqIO

def rmPE(read1,read2,adaptor1,adaptor2,min_length):
res_1 = rmSE(read1,adaptor1,min_length)
res_2 = rmSE(read2,adaptor2,min_length)
if res_1 and res_2:
return res_1,res_2
else:
return False

def rmSE(read,adaptor,min_length):
seq = read.seq
seed_len = 6
a_len = len(adaptor)
seq_len = len(seq)
for i in range(a_len - seed_len):
seed = adaptor[i:i+seed_len]
pos = 0
while(pos < seq_len):
find_pos = seq.find(seed,pos)
if find_pos > 0:
mistaken_count = 0
_b = find_pos
_e = find_pos + seed_len
while(_b >= 0 and i >= find_pos - _b):
if adaptor[i - find_pos + _b] != seq[_b]:
mistaken_count += 1
if mistaken_count > 3:
break
_b -= 1
else :
while(_e < seq_len and i - find_pos + _e < a_len):
if adaptor[ i - find_pos + _e ] != seq[_e]:
mistaken_count += 1
if mistaken_count > 3:
break
_e += 1
else:
if find_pos - i > min_length:
return  read[:find_pos-i]
else :
return False
pos = find_pos + 1
else:
break
return read

def rmAdaptor(argv):
argv.pop(0)
type = argv.pop(0)
if type=='PE':
read1_file,read2_file,adaptor1,adaptor2,out_prefix,min_length = argv
read2_records = SeqIO.parse(open(read2_file),'fastq')
read1_out = open( '%s.1.fq'%out_prefix,'w' )
read2_out = open( '%s.2.fq'%out_prefix,'w' )
for read1 in SeqIO.parse(open(read1_file),'fastq'):
read2 = read2_records.next()
rmPE_res = rmPE(read1,read2,adaptor1,adaptor2,min_length)
if rmPE_res:
read1_out.write(rmPE_res[0].format('fastq'))
read2_out.write(rmPE_res[1].format('fastq'))
elif type=='SE':
reads_file,adaptor,out_prefix,min_length = argv
reads_out = open( '%s.single.fq'%out_prefix,'w' )
for reads in SeqIO.parse(open(reads_file),'fastq'):
rmSE_res = False
if re.search('[\s\/](\d)',reads.id).group(1) == '1':
rmSE_res = rmSE(reads,adaptor1,min_length)
elif re.search('[\s\/](\d)',reads.id).group(1) == '2':
rmSE_res = rmSE(reads,adaptor2,min_length)
if rmSE_res:
reads_out.write(rmSE_res.format('fastq'))

if __name__ == '__main__':
rmAdaptor(sys.argv)