Оптимизация моего скрипта, который ищет большой сжатый файл - PullRequest
Новая
       48

Оптимизация моего скрипта, который ищет большой сжатый файл

2 голосов
/ 21 февраля 2020

Я снова здесь! Я хотел бы оптимизировать мой bash скрипт, чтобы уменьшить время, затрачиваемое на каждый l oop. По сути, это:

  • получение информации из tsv
  • с использованием этой информации для поиска с помощью awk в файле
  • , печати строки и ее экспорта

Мои проблемы: 1) файлы представляют собой сжатые файлы размером 60 ГБ: мне нужно программное обеспечение для его распаковки (сейчас я пытаюсь распаковать его, не уверен, что у меня будет достаточно места) 2) В любом случае, долго разбираться в этом

Мои идеи улучшить его:

  1. 0) как я уже сказал, если возможно, я распакую файл

  2. используя GNU параллельно с parallel -j 0 ./extract_awk_reads_in_bam.sh ::: reads_id_and_pos.tsv, но я не уверен, что он работает как положено? Я сокращаю время на исследование с 36 минут до 16, так что всего в 2,5 раза? (У меня 16 ядер)

  3. Я думал (но это может быть избыточно с GNU?), Чтобы разделить мой список информации, чтобы просмотреть несколько файлов, чтобы запустить их параллельно

  4. сортировка файла bam по имени чтения и выход из awk после того, как найдено 2 совпадения (не может быть больше 2)

Вот остальные мои bash Сценарий, я действительно открыт для идей по его улучшению, но я не уверен, что я суперзвезда в программировании, так что, возможно, простота поможет? :)

Мой bash скрипт: 

#/!bin/bash
while IFS=$'\t' read -r READ_ID_WH POS_HOTSPOT; do
echo "$(date -Iseconds) read id is : ${READ_ID_WH} with position ${POS_HOTSPOT}" >> /data/bismark2/reads_done_so_far.txt
echo "$(date -Iseconds) read id is : ${READ_ID_WH} with position ${POS_HOTSPOT}"
samtools view -@ 2 /data/bismark2/aligned_on_nDNA/bamfile.bam | awk -v read_id="$READ_ID_WH" -v pos_hotspot="$POS_HOTSPOT" '$1==read_id {printf $0 "\t%s\twh_genome",pos_hotspot}'| head -2 >> /data/bismark2/export_reads_mapped.tsv
done <"$1"

Мой tsv-файл имеет такой формат: 

READ_ABCDEF\t1200

Большое спасибо ++

1 Ответ

1 голос
/ 21 февраля 2020

TL; DR

Ваш новый скрипт будет: 

#!/bin/bash
samtools view -@ 2 /data/bismark2/aligned_on_nDNA/bamfile.bam | awk -v st="$1" 'BEGIN {OFS="\t"; while (getline < st) {st_array[$1]=$2}} {if ($1 in st_array) {print $0, st_array[$1], "wh_genome"}}'

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

samtools view -@ 2 "$bam" | grep -f <(awk -F$'\t' '{print $1}' "$1") > "$sam"

Здесь вы получаете все чтения с samtools и ищете все термины, которые встречаются в -f параметр grep. Этот параметр является файлом, который содержит первый столбец входного файла поиска. В результате получается файл sam, содержащий только те чтения, которые перечислены в файлах ввода поиска. 

awk -v st="$1" 'BEGIN {OFS="\t"; while (getline < st) {st_array[$1]=$2}} {print $0, st_array[$1], "wh_genome"}' "$sam"

Наконец, используйте awk для добавления дополнительной информации:

  1. Открыть входной файл поиска с awk в начале и считайте его содержимое в массив (st_array)
  2. Установите разделитель поля вывода на табулятор
  3. Перейдите в файл sam и добавьте дополнительную информацию из предварительно заполненного массива.

Я предлагаю эту схему, потому что я чувствую, что grep быстрее, чем awk для выполнения поиска, но тот же результат может быть получен только с помощью awk : 

samtools view -@ 2 "$bam" | awk -v st="$1" 'BEGIN {OFS="\t"; while (getline < st) {st_array[$1]=$2}} {if ($1 in st_array) {print $0, st_array[$1], "wh_genome"}}'

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

В любом случае вам необходимо перечитать файл более одного раза или распакуйте его перед началом работы с ним.

...