BioPython: обработка сырых RNAseq-считываний (фильтрация качества и обрезка)

Question

У меня есть необработанный файл fastq.gz, который я пытаюсь предварительно обработать с помощью Biopython перед выравниванием. В конечном итоге я хотел бы удалить показания низкого качества, обрезать polyA-хвосты, обрезать адаптеры с использованием нечеткого соответствия и, наконец, удалить считывания, которые не удовлетворяют требованию длины после всей указанной предварительной обработки. Также было бы аккуратно указать, сколько чтений удовлетворяет критериям фильтрации на каждом шаге. Я играл с этим сценарием на Biopython, но имел небольшой успех.

# coding: utf-8 

''' 
process_reads.py 
March 31, 2016 

Convert FQ files to unaligned and tagged BAMs 
I: compressed FQ 
O: compressed, filtered FQ 
''' 

import gzip, statistics 
from Bio import SeqIO, pairwise2 

def get_stats(reads): 
    sizes = [len(r) for r in reads] 
    print("Total reads: %i" % len(sizes)) 
    print("Mean read length: %i" % statistics.mean(sizes)) 
    print("Max. read length: %i" % max(sizes)) 
    print("Min. read length: %i" % min(sizes)) 

def quality_filter(reads, qual): 
    return (r for r in reads if min(r.letter_annotations["phred_quality"]) >= qual) 

def trim_polyA(records, numA, minLen): 
    for record in records: 
     if len(record) < minLen: continue 
     record = record.seq.split("A"*numA, 1)[0] 
     yield record 

def _remove_adaptor(seq, region, right_side=True): 
    if right_side: 
     try: 
      pos = seq.find(region) 
     # handle Biopython SeqRecords 
     except AttributeError: 
      pos = seq.seq.find(region) 
     return seq[:pos] 
    else: 
     try: 
      pos = seq.rfind(region) 
     # handle Biopython SeqRecords 
     except AttributeError: 
      pos = seq.seq.rfind(region) 
     return seq[pos+len(region):] 

def trim_adaptor(seq, adaptor, num_errors, right_side=True): 
    gap_char = '-' 
    exact_pos = str(seq).find(adaptor) 
    if exact_pos >= 0: 
     seq_region = str(seq[exact_pos:exact_pos+len(adaptor)]) 
     adapt_region = adaptor 
    else: 
     seq_a, adaptor_a, score, start, end = pairwise2.align.localms(str(seq), 
                     str(adaptor), 
                     5.0, -4.0, -9.0, -0.5, 
                     one_alignment_only=True, 
                     gap_char=gap_char)[0] 
     adapt_region = adaptor_a[start:end] 
     seq_region = seq_a[start:end] 
    matches = sum((1 if s == adapt_region[i] else 0) for i, s in enumerate(seq_region)) 
    # too many errors -- no trimming 
    if (len(adaptor) - matches) > num_errors: 
     return seq 
    # remove the adaptor sequence and return the result 
    else: 
     return _remove_adaptor(seq, seq_region.replace(gap_char, ""), 
       right_side) 

def process_reads(fq, qual, adapt, numA, minLen): 
    with gzip.open(fq) as f: 
     rawReads = SeqIO.parse(f, "fastq") 
     # get_stats(rawReads) # When I run this, everything downstream fails.. 
     qualFil = quality_filter(rawReads, qual) # I think this work fine. 
     trimmedPoly = trim_polyA(qualFil, numA, minLen) 
     trimmedAdap = trim_adaptor(trimmedPoly, adapt, 2) 

     # count = SeqIO.write(trimmedAdap, "good_quality.fastq", "fastq") 
     # print(count) 

# TEST PROCESSING 
fq = "test/TAAGGCGA_2.fq.gz" 
process_reads(fq, qual=50, adapt="AAGCAGTGGTATCAACGCAGAGTGAATGGG", numA=6, minLen=20) 

Я считаю, что фильтр качества и обрезка polyA работают правильно, но я не могу заставить адаптеры отрезать. Я также написал функцию под названием get_stats, которая должна возвращать среднюю длину и общее число просмотров. Буду признателен за любую помощь!

Answer 1

Мне нравится, так как right_side всегда установлен на True. Таким образом, ваш метод _remove_adaptor только попытаться извлечь адаптер с правой стороны:

 
... 
if right_side: 
     try: 
      pos = seq.find(region) 
     # handle Biopython SeqRecords 
     except AttributeError: 
      pos = seq.seq.find(region) 
     return seq[:pos] 

Но если отрезать полиА с правой стороны (удерживая левую часть последовательности), вы все еще есть последовательность адаптера справа? Я предполагаю, что вы хотите расщепить адаптер с левой стороны.

Некоторые примеры чтения были бы полезными.