Соответствующие части строк в файле (python)

Question

В настоящее время у меня есть список генов в файле. В каждой строке есть хромосома со своей информацией. Такая запись выглядит как:

NM_198212 chr7 + 115926679 115935830 115927071 11593344 2 115926679, '115933260', 115927221, '115935830',

Последовательность для хромосомы начинается у основания 115926679 и продолжается до (но не включая) основания 115935830

Если нам нужна последовательность сращивания, мы используем экзоны. Первый простирается от 115926679 до 155927221, а второй - от 115933260 до 115935830.

Тем не менее, я столкнулся с проблемой при дополнительной последовательности, такой как:

NM_001005286 chr1 - 245941755 245942680 245941755 245942680 1 245941755, '245942680'

Поскольку столбец 3 является "-", эти координаты относятся к антисмысловой цепи (дополнение к цепи). Первая база (жирным шрифтом) соответствует последней базе смысловой цепи (курсивом). Поскольку файл имеет только смысловую позицию, мне нужно попытаться перевести координаты антисмысловой цепи в смысловую цепь, выбрать правильную последовательность и затем дополнить ее в обратном порядке.

Тем не менее, я программирую только около полугода и не знаю, как начать делать это.

Я написал регулярное выражение:

'(NM_\d+)\s+(chr\d+)([(\+)|(-)])\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+),(\d+),s+(\d+),(\d+),'

но теперь я не уверен, как запустить эту функцию... Если кто-нибудь может помочь мне вообще начать с этим, возможно, заставит меня понять, как это сделать, я был бы очень признателен.

ОК: предположим, что это хромосома 25:

AAAAAAAAAACCCCCCCCCCTTTTTTTTTTGGGGGGGGGG

(есть 10 каждого персонажа).

Теперь: если я ищу не сплайсированный ген на: chr25 + 10 20

Затем ген начинает с позиции 10 (начиная с 0) и поднимается до позиции 20, но не включает ее.

CCCCCCCCCC

Это просто. Это очень хорошо соответствует нарезке строк Python.

Это более запутанно, если я дам вам:

chr25 - 10 20

То, что у вас есть, это положительная нить. Но этот ген находится на отрицательной (комплементарной) цепи. Помните, как хромосома выглядит как двухцепочечная:

AAAAAAAAAACCCCCCCCCCTTTTTTTTTTGGGGGGGGGG
TTTTTTTTTTGGGGGGGGGGAAAAAAAAAACCCCCCCCCC

Мы ищем ген на нижней цепи. То есть мы считаем от 0, начиная справа. Пронумеруйте верхнюю нить слева, а нижнюю нить справа. Итак, что я хочу здесь, это AAAAAAAAAA.

Подвох в том, что я даю вам только верхнюю прядь. Я не дам тебе нижнюю нить. (Вы можете сгенерировать себя из верхней нити - но учитывая, насколько она велика, я советую против этого.)

Так что вам нужно конвертировать координаты. На нижней нити основание 0 (крайняя справа C) противоположно основанию 39 на верхней нити. База 1 против базы 38. База 2 против случая 37. (Важный момент: обратите внимание, что происходит, когда вы добавляете эти два числа вверх - каждый раз.) Таким образом, база 10 против базы 29, а база 19 против базы 20.

Итак: если я хочу найти базу 10-20 на нижней нити, я могу посмотреть на базу 20-29 на верхней (а затем дополнить ее обратно).

Мне нужно выяснить, как перевести координаты на нижней нити в эквивалентные координаты на нижней нити. Да: это очень запутанно

Я попробовал оригинальный ответ Вероники:

fields = line.split(' \t')
geneID, chr, strand = fields[:2]
start = int(fields[3])
end = int(fields[4])
if strand == '-':
    start,end = -(start + 1), -(end + 1) # this part was changed from original answer.

который находится на правильном пути, но этого недостаточно. Это займет 10 и 20, и превратить его в 20 и 10.

И я знаю, что могу изменить строку, выполнив это:

r = s[::-1]
bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
l = list(r)
o = [bc[base] for base in l]
     return ''.join(o)

Под редакцией! Это выглядит правильно?!

fp2 = open('chr22.fa', 'r')
fp = open('chr22.fa', 'r')
for line in fp2:
    newstring = ''
    z = line.strip()
    newstring += z
for line in fp:
    fields = line.split('\t')
    gene_ID, chr, strand = fields[:2]
    start = int(fields[3])
    end = int(fields[4])
    bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'a': 't', 't': 'a', 'c':'g', 'g':'c', 'N':'N', 'n':'n'}
    l = list(newstring)        
    if strand == '+':
        geneseq = ''.join([bc[base] for base in l[start:end]]) 
    if strand == '-':
        newstart, newend = -(start + 1), -(end + 1)
        genseq = ''.join([bc[base] for base in l[newstart:newend]]) 

Answer 1

Если вы нарезаете строку с отрицательным числом, Python считает в обратном направлении от конца.

complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
s = "AAAAAAAAAACCCCCCCCCCTTTTTTTTTTGGGGGGGGGG"
"".join(complement[c] for c in s[-20:-10])

---

РЕДАКТИРОВАТЬ:

Ваше редактирование выглядит как раз для меня, да. Я очень плохо проверяю ошибки на заборе, но вы все равно лучше, чем я!

Я переписал твой код, чтобы он был более Pythonic, но не изменил ничего существенного.

bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N':'N'}

f = open('chr22.fa')
l = ''.join(line.strip() for line in f)
f.seek(0)

for line in f:
    fields = line.split('\t')
    gene_ID, chr, strand = fields[:2]

    start = int(fields[3])
    end = int(fields[4])

    if strand == '-':
        start, end = -(start + 1), -(end + 1)

    geneseq = ''.join(bc[base.upper()] for base in l[start:end])

Answer 2

Как отмечено в моем комментарии к вопросу, это выглядит как очень странный формат файла, таким образом, мое первоначальное замешательство.

Примечание: если это один из стандартных форматов биологических файлов, вам, вероятно, лучше использовать Biopython или что-то подобное для его анализа.

Если вы хотите выполнить собственный синтаксический анализ, регулярные выражения по-прежнему кажутся неправильным путем - их трудно читать и не нужно с простым файлом, разделенным пробелами и табуляцией. Я предполагаю, что вы проанализировали файл фаста хромосом и у вас есть последовательности всех хромосом как chrom_sequences имя: seq словарь, а также что у вас есть reverse_complement функция (обе эти вещи легко реализуются вручную, но, вероятно, лучше с биопионом).

fields = line.split(' ')  # or '\t' instead of ' ' if the file is tab-separated
gene_ID, chr, strand = fields[:2]
start, end = [int(x) for x in fields[3:5])
this_chromosome_seq = chrom_sequences[chr]
# if strand is +, just take the sequence based on the start-end position
if strand == '+':
  # be careful about off-by-one errors here - some formats give you a 1-based position, 
  #  other formats make it 0-based, and they can also be end-inclusive or end-exclusive
  gene_sequence = this_chromosome_seq[start:end]
# if your coordinates really are given as antisense strand coordinates when strand is -,
#  you just need to subtract them from the chromosome length to get sense-strand coordinates,
#  (switching start and end so they're still in smaller-to-larger order),
#  and then reverse-complement the resulting sequence. 
if strand == '-':
  chrom_length = len(this_chromosome_seq)
  # again, be careful about off-by-one errors here!
  new_start,new_end = chrom_length-end, chrom_length-start
  gene_sequence = reverse_complement(this_chromosome_seq[new_start:new_end])

Оригинальный ответ, на самом деле не делал то, что просили:

Если вы просто хотите получить начальную / конечную позиции, сделайте что-то вроде этого:

fields = line.split(' ')  # or '\t' instead of ' ' if the file is tab-separated
gene_ID, chr, strand = fields[:2]
start = int(fields[3])
end = int(fields[4])
if strand == '-':
  start,end = end,start

Затем вам нужно будет проанализировать ваш файл fasta, чтобы получить последовательность этих начальных и конечных координат, и дополнить ее, если strand=='-', Опять же, я думаю, что Biopython может сделать большую часть этого для вас.

Еще одно замечание - Biostar - это хороший сайт типа StackExchange, специально предназначенный для биоинформатики, ответы на который вы можете получить там.

Answer 3

Я не понял проблемы домена, но похоже, что вы пытаетесь втиснуть слишком много в одно регулярное выражение. Попробуйте разбить его на более простые подзадачи, например, следующие (в псевдокоде):

if third column is a '+'
    parseRegularSequence()
else
    parseComplementarySequence()

Answer 4

Я полагаю (особенно потому, что ваши файлы большие), что было бы НАМНОГО проще читать и писать прямо из файлового буфера.

Допустим, вы уже проанализировали свой заголовочный файл. Разобранная вами строка выглядит так:

line = "NM_001005286 chr1 - 245941755 245942680 245941755 245942680 1"

Затем вы определяете начальную / конечную позиции (в антисмысловых координатах):

name, chromosome, direction, start, end = line[:5]

Затем сделайте следующее:

#Open up the file `chr1.txt`.
f = open("chr1.txt", "r")

#Determine the read length.
read_len = end - start

#Seek to the appropriate position (notice the second argument, 2 -- this means
#seek from the end of the file)
f.seek(-end, 2)

#Read the data
sense_seq = f.read(read_len)

После этого нужно просто преобразовать последовательность в антисмысловую.

Легко:)