linux - BASH 使用从另一个文件中的一列传输的值递归地制作堆积文件

Question

我正在尝试使用 samtools 从两个文件 File1 和 File2 制作 pileup 文件。

我按染色体拆分了 File1 和 File2，导致有 44 个文件按照以下格式命名:

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

其中 ${c} 是 1 到 22 之间的数字，$TISSUE 是结肠或肌肉——22 条染色体用于结肠，22 条用于肌肉。 IE。; 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 的新手，我基本上没有接受过计算机科学培训。

Answer 1

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。循环。

我不明白你一开始是如何拆分这些文件的，或者什么是 pileup，但我怀疑如果你直接处理 File1，这可以大大简化。和 File2相反。

您也可以避免外循环，只对所有 *_ONLY 运行 Awk直接存档。您可以从 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]从文件名中获取染色体标识符部分。