Я хотел бы разделить файл sam на несколько файлов sam в соответствии с информацией о штрих-коде.И информация о штрих-коде запроса приведена в другом файле.
$ cat barcode.list
ATGCATGC
TTTTAAAA
GGGGCCCC
CGCGATGA
AAGGTTCC
....
Простой скрипт bash, приведенный ниже, может достичь цели.
barcode_list=./A_barcode.csv
input_bam=./A_input.bam
splited_dir=./splited_sam/A
filtered_dir="./filterd_sam/A"
mkdir -p ${splited_dir} ${splited_dir}
header=$(samtools view -H ${input_bam})
samtools view {input.bam} | LC_ALL=C fgrep -f <(cat ${barcode_list}) | awk -v header="${header}" -v outdir="${splited_dir}" '{barcode=substr($0,index($0, "\tCB:Z:")+6,18);if (!header_printed[barcode]++) {print $0 >> outdir"/"barcode".sam"}}'
for sam in ${output_dir};do samtools view -q 30 -m 1 ${sam} -O bam -o ${filtered_dir}/$(basename ${sam} "sam")"bam";done
Примечание: В новый файл будут записаны только штрих-коды, которые есть как в файле barcode_list, так и в файле input_bam.
Но я хочупереписать скрипт в конвейер sankemake для лучшего масштабирования.Решение, которое я попробовал, показано ниже.
Я не знаю, как назначить имя входного файла на последнем шаге всех правил , rule all в этом примере.Потому что они определяются как input_bam, так и input_barcode файлом.Между тем, без знания splited_sam имени файла, я не могу перейти к следующему шагу.
SAMPLES = ["A", "B", "C", "D"]
# BARCODE = ???
rule all:
input:
splited_sam_dir = expand("splited_sam/{sample}", sample=SAMPLES)
rule split_sam:
input:
bar = "{sample}_barcode.csv",
bam = "{sample}_input.bam"
output:
splited_sam_dir = "splited_sam/{sample}"
shell:
"""
header=$(samtools view -H {input.bam})
samtools view {input.bam} | LC_ALL=C fgrep -f <(cat {input.bar}) | awk -v header="$header" -v outdir="{output.splited_sam_dir}" '{{barcode=substr($0,index($0, "\tCB:Z:")+6,18);if (!header_printed[barcode]++) {{print $0 >> outdir"/"barcode".sam"}}}}
"""
rule filter_sam:
# ??? don't know the input file name...