0

我正在尝试使用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 条染色体代表肌肉。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
.
.
.

这些文件由两列组成，第一列只显示染色体编号，第二列是该染色体上的位置。IE;

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、 wherea=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 的新手，而且我基本上没有受过计算机科学培训。

4

2 回答 2

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]以从文件名中获取染色体标识符部分。

于 2017-12-31T22:22:30.957 回答
0

这最终为我工作：

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

谢谢您的帮助！

于 2018-01-12T15:27:22.947 回答