Как найти файл для всех позиций в диапазоне в Perl или Python?

Question

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

gene_id start stop 
grf4      1245  1365
fgt89     3089  4524
tig3      45600 46800

Другойсодержит позиции и информацию об этих позициях, например:

position   nucleotide  %support
3980         T          98%
456          C          78%
45900        G          100%
4234         C          70%

Я хотел бы сгенерировать файл, в котором есть строка для каждого gene_id и вся информация из второго файла, который попадает в этот ген, например:

gene_id  start  stop position nucleotide support
fgt89    3089    4524  3980    T         98%
fgt89    3089    4524  4234    C         70%
tig3     45600   46800 45900   G         100% 

Я не уверен, что лучший способ сделать это, любые предложения приветствуются!Файлы довольно большие, поэтому скорость - проблема.Спасибо!

Answer 1

Изменить текст в данные

Сначала я бы предложил обработать ваши файлы и сохранить их как базу данных для облегчения дальнейшей обработки. Вы можете извлечь данные из этих файлов, например. следующими способами:

file1 = open('some_filename.txt')
file1_lines = [filter(len, line.split(' ')) for line in file1.readlines()]

и file1_lines будет содержать что-то похожее на это:

[['gene_id', 'start', 'stop'], ['grf4', '1245', '1365'], ['fgt89', '3089', '4524'], ['tig3' , '45600', '46800']]

, что означает, что у вас есть данные, которые вы можете вставить в базу данных.

Вставить данные в базу данных

Теперь вы должны вставить данные в базу данных, чтобы их можно было легко обработать с помощью соответствующих инструментов (в данном случае СУБД). Я предлагаю использовать SQLite (см. Встроенный модуль sqlite3 для Python ). Если вы не знаете SQL, задайте ему конкретные вопросы. Но я считаю, что база данных - это форма, наиболее подходящая для вашей задачи.

Извлечение необходимых данных из базы данных

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

SELECT `ranges`.`gene_id`, `ranges`.`start`, `ranges`.`stop`,
    `nucleotides`.`position`, `nucleotides`.`nucleotide`, `nucleotides`.`support`
FROM `ranges`
JOIN `nucleotides`
ON (
    (`ranges`.`start` < `nucleotides`.`position`)
    AND
    (`nucleotides`.`position` < `ranges`.`start`)
)
ORDER BY `ranges`.`gene_id` ASC, `nucleotides`.`position` ASC

Что даст вам нужные данные в порядке gene_id по возрастанию, затем position по возрастанию.

Единственная задача, которая останется, - это вывод данных в файл. Для этого я предлагаю использовать format().

Основная информация

Итак, ваша задача состоит в основном из следующих частей:

1. Превратить текст (содержимое файлов) в данные (списки переменных).
2. Создать базу данных с правильной структурой.
3. Вставка данных из файлов в базу данных (убедитесь, что целые числа и числа с плавающей запятой хранятся как таковые, а не как строки).
4. Сделайте запрос к базе данных, как я упоминал выше.
5. При необходимости вывести данные в файл в виде текста (используя форматирование строки - format()).

Это помогло? Есть вопросы?

Answer 2

База данных вне курса - это идеальное решение, опубликованное Tadeck, но вот что вы можете попробовать

>>> def Test():
    fin1=open("file1.txt") #File as per your First Table
    fin2=open("file2.txt") #File as per your Second Table
    fin1.readline()        #Skip the Header
    fin2.readline()        #Skip the Header
    #Sort The First list and create an Iterator
    data1=iter(sorted([[f.split()[1],f.split()[2],f.split()[0]] 
               for f in fin1.xreadlines()], key=operator.itemgetter(0)))  
    #Sort The Second List and create an Iterator
    data2=iter(sorted([f.split() for f in fin2.xreadlines()],
               key=operator.itemgetter(0))) 
    #Print The Header
    print "{0:10}{1:10}{2:10}{3:10}{4:10}{5:10}".format("gene_id",
                                                        "start",
                                                        "stop",
                                                        "position",
                                                        "nucleotide",
                                                        "support") 
    try:
        v1=data1.next() #Read First Item from First List
        v2=data2.next() #Read Second Item from First List
        while True: #Loop Until One of the List has reached the end
            #If the Position is greater than stop range (from first list), 
            #read the next item from the first list
            if v2[0] > v1[1]: 
                 v1=data1.next() 
            #If the Position is greater or equal than the start range 
            #(We are in the range)
            elif v2[0] >= v1[0]: 
                #Format and Print it
                print "{0:10}{1:10}{2:10}{3:10}{4:10}{5:10}".format(v1[2],
                                                                    v1[0],
                                                                    v1[1],
                                                                    v2[0],
                                                                    v2[1],
                                                                    v2[2]) 
                #Read the Next Item From the Second List
                v2=data2.next() 
            #Not in any Range so Read the Next Item From the Second List
            else: v2=data2.next() 
    except StopIteration:
        None


>>> Test()
gene_id   start     stop      position  nucleotidesupport   
fgt89     3089      4524      3980      T         98%       
fgt89     3089      4524      4234      C         70%       
tig3      45600     46800     45900     G         100%      
>>> 

Answer 3

Простое решение Perl может выглядеть так:

#!/usr/bin/perl
use warnings;
use strict;

open my $F2, '<', '8410347-file2.txt' or die $!;
my %at_position;
while (<$F2>) {
    my ($position, $nucleotide, $support) = split;
    warn "Duplicate position $position\n" if exists $at_position{$position};
    $at_position{$position} = {nucleotide => $nucleotide,
                               support    => $support};
}
close $F2;

open my $F1, '<', '8410347-file1.txt' or die $!;
while (<$F1>) {
    my ($id, $start, $stop) = split;
    my @matches = grep $_ >= $start && $_ <= $stop, keys %at_position;
    for my $match (@matches) {
        print join("\t",
                   $id,
                   $start,
                   $stop,
                   $match,
                   $at_position{$match}{nucleotide},
                   $at_position{$match}{support},
                  ), "\n";
    }
}

Если ваши файлы большие, этот подход будет слишком медленным, и потребуется некоторая более сложная структура для захвата интервалов (например, см. Календарь Ленивого Менеджера ).