Как сделать поиск мотивов с помощью Python?

Question

Я пытаюсь проверить мотив привязки nrf2, используя регулярное выражение с python.Я сделал это с R, используя JASPAR2018 PWM, но из-за нескольких проблем с JASPAR.

Я хочу переделать его, используя python.

Попытка

from Bio import SeqIO
from itertools import islice
import pandas as pd

#Creating Reverese Complements
def reverseComp(Seq):
    seq = Seq.upper()
    d = {'A':'T', 'T':'A', 'G':'C', 'C':'G'}
    try:
        seq = seq[::-1]
        rc_seq = "".join([d[nuc] for nuc in seq])
    except KeyError:
        return "Not Viable DNA Seq"
    return rc_seq

def genSeq(genome_path, chrom, chromstart, chromend):
    if bool(re.search('gz', genome_path)) | bool(re.search('fa', genome_path)) | bool(re.search('fasta', genome_path)):
        if bool(re.search('gz', genome_path)) == True:
            genome = SeqIO.parse(gzip.open(genome_path, 'rt'),'fasta')
            identifiers = [seq_record.id for seq_record in genome]
            seq_gen = next(islice(genome, identifiers.index(chrom) , None))
            seq = str(seq_gen.seq[chromstart:chromend])
        else:
            genome = SeqIO.parse(open(genome_path),'fasta')
            identifiers = [seq_record.id for seq_record in genome]
            seq_gen = next(islice(genome, identifiers.index(chrom)+1 , None))
            seq = str(seq_gen.seq[chromstart:chromend])
    elif bool(re.search('2bit', genome_path)):
        tbGenome = tbr.TwoBitFile(genome_path)
        seq = tbGenome[chrom][chromstart:chromend]
    else:
        raise Exception('File type not recognized')
    return (seq).upper()

pat = "[AGC]TGA[CTG][ATCG][CAT][AGT]GC[ATCG]"
pattern = re.compile(pat)

motifDF = []
motifQuant = []
with open('/Users/kalyanidhusia/Desktop/nrf2_R/ENCFF126HBJ.bed') as f:
    for line in f:
        peak = list(line.split())
        seq = genSeq('hg19.fa', peak[0], int(peak[1]), int(peak[2]))
        rSeq = reverseComp(seq)
        sequences = []
        for result in re.finditer(pattern, seq):
            sequences.append("".join(result.groups()))
        for result in re.finditer(pattern, rSeq):
            sequences.append("".join(result.groups()))
        if len(sequences) > 0:
            seqs = pd.DataFrame({'binding':sequences, 'chrom':peak[0], 'chromstart':peak[1], 'chromend':peak[2]})
            motifDF.append(seqs)
            motifQuant.append([peak[0], peak[1], peak[2], len(seqs), len(seq)])
search_reg = pd.concat(motifDF)
names = ['chrom', 'chromstart', 'chromend', 'numOfMatches', 'lenSeq']
dist_reg = pd.DataFrame(motifQuant, columns=names)

Ошибка

Это ошибка, которую я получаю:

ipython-input-3-2e7ebdf92205> в genSeq(genome_path, chrom, chromstart, chromend) 25 идентификаторов = [seq_record.id для seq_record в геноме] ---> 26 seq_gen = next (islice (геном, identifiers.index (chrom) +1, нет)) 27 seq = str(seq_gen.seq [chromstart: chromend]) 28 elif bool (re.search ('2bit', genome_path)): StopItered:

Как мне решить эту проблему?

Answer 1

К вышеупомянутой проблеме я смог ее решить, немного подправив свой код. Вот решенный пример для вас, ребята, и моя проблема с кодом ниже:

motif = '[REGULAR_EXPRESSION_FOR_YOUR_MOTIF]'
regBS = re.compile(motif)
motifDF = []
motifQuant = []
genome = tbr.TwoBitFile('/Path_to_your_genomefile_in_2bit.2bit/')
with open('/Path_to_your.bedfile/') as f:
    for line in f:
        if line.startswith('track') == False:
            peak = list(line.split())
            seq = (genome[peak[0]][int(peak[1]):int(peak[2])]).upper()
            rSeq = reverseComp(seq)
            sequences = []
            sequences.extend(re.findall(regBS, seq))
            sequences.extend(re.findall(regBS, rSeq))
            if len(sequences) > 0:
            seqs = pd.DataFrame({'binding':sequences, 'chrom':peak[0],'chromstart':peak[1], 'chromend':peak[2], 'NR':'NRF2'})
                motifDF.append(seqs)
                motifQuant.append([peak[0], peak[1], peak[2], len(seqs), len(seq)])
         search_reg = pd.concat(motifDF)
names = ['chrom', 'chromstart', 'chromend', 'numOfMatches', 'lenSeq']
dist_reg = pd.DataFrame(motifQuant, columns=names)
dist_reg.head()
n = 5
x = [len(i[6+n:-6-n]) for i in search_reg['binding']]

Этот код генерирует пиковые последовательности, которые мне нужны, и сохраняет их в search_reg [binding], но он также хранит вместе с ним разделенное пробелами число. Мне нужно хранить их в двух разных столбцах. Есть предложения?