BASH 使用从另一列中的一列通过管道传输的值递归创建堆积文件

Question

我正在尝试使用两个文件File1和File2中的samtools制作堆积文件。

我按染色体拆分了 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")，我需要从第 2 列中获取位置，将其称为"x"，并使用它来获取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循环。

我不明白你首先是如何拆分这些文件的，或者堆积是什么，但我怀疑如果你只是直接处理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]，以从文件名中获取染色体标识符部分。

Answer 2

这就是最终对我有用的：

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

感谢您的帮助！