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

Я пытаюсь создать файлы наложения, используя samtools из двух файлов, File1 и File2.

Я разделил File1 и File2 по хромосомам, в результате чего 44 файла были названы в следующем формате:

chr${c}.${TISSUE}_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY

где ${c} — число от 1 до 22, а $TISSUE — либо толстая кишка, либо мышца: 22 хромосомы для толстой кишки и 22 для мышц. то есть; chr1.colon_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY chr2.colon_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY

.
.
.

chr22.colon_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY
chr1.muscle_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY
.
.
.

Эти файлы состоят из двух столбцов, первый просто показывает номер хромосомы, а второй столбец — положение на этой хромосоме. то есть;

head chr2.colon_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY 
chr2 103977
chr2 112051
chr2 126199
chr2 146288
chr2 147797
chr2 147822
chr2 148548
chr2 148525
chr2 158189
chr2 158188

Для каждой строки в файле (например, "chr2.colon_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY") мне нужно взять позицию, назвать ее «x» из столбца 2 и использовать ее для получения диапазона a-b, где a=x-5 и b=x+5. Затем я подключу эти значения к следующему скрипту:

samtools mpileup -f [REFERENCE GENOME] File1 File2 -r chr${c}:a-b

Например, предположим, что я смотрю на хромосому 2, позиция 103977 (строка 1 выше). Тогда мой сценарий будет

samtools mpileup -f [REFERENCE GENOME] File1 File2 -r chr2:103972-103982

В общем, это петля внутри петли внутри петли. Что-то типа,

for t in $(colon, muscle)
do
  for c in $seq (1 22)
  do
    for item (or maybe row?) in 
      chr${c}.${t}_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY
    do
      awk '{print $2}' | something something something 
      x= position in col 2, a=x-5 b=x+5
      samtools mpileup -f [REFERENCE GENOME] File1 File2 -r chr${c}:a-b
    done 
  done
done
...

Заранее спасибо. Я новичок в работе с Linux и практически не имею никакого образования в области компьютерных наук.


bash linux loops awk samtools
person Emm Gee    schedule 31.12.2017    source источник
comment
Привет, пожалуйста, используйте редактирование и используйте функцию кода (фигурные скобки) в редакторе сообщений для удобочитаемости. Вопрос не читается, пожалуйста, организуйте его. Я рекомендую вам посмотреть на другие вопросы, чтобы узнать, как написать правильный вопрос. ГЛ :)   -  person Blacky    schedule 31.12.2017

Ответы (2)


arrow_upward
1
arrow_downward

Awk обрабатывает строку за раз, поэтому я бы выбрал что-то вроде

for t in colon muscle; do
    for c in $(seq 1 22); do
        awk '{ print $2-5 "-" $2+5 }' chr${c}.${t}_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY |
        while read -r range; do
            samtools mpileup -f [REFERENCE GENOME] File1 File2 -r chr${c}:$range
        done 
    done
done

Другими словами, Awk обрабатывает файл целиком и передает по одной строке вывода в последний цикл while read -r range.

Я вообще не понимаю, как вы разделяете эти файлы или что такое пайлап, но я подозреваю, что это можно было бы значительно упростить, если бы вы вместо этого работали непосредственно с File1 и File2.

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

awk '{ print $1 ":" $2-5 "-" $2+5 }' *_ONLY |
while read -r chrrange; do
    samtools mpileup -f [REFERENCE GENOME] File1 File2 -r "$chrrange"
done

Если вы не можете использовать $1 напрямую, попробуйте split(FILENAME, f, /\./) и напечатайте f[1], чтобы получить часть идентификатора хромосомы из имени файла.

person tripleee    schedule 31.12.2017

arrow_upward
0
arrow_downward

Вот что в итоге сработало для меня:

module load SAMtools

awk '{print $1, $2-5 "-" $2+5}' FILE PATH |\
while read chrom range
do

    samtools mpileup -f /REFERENCE GENOME\
            /${chrom}.COLON BAM FILE\
            /${chrom}.MUSCLE BAM FILE\
            -r $chrom:$range -o ${chrom}.colon.${range}.pileup

done

Спасибо за помощь!

person Emm Gee    schedule 12.01.2018