我有几个后缀为1.ht2l至.8.ht2l的基因组文件
bob.1.ht2l
...
bob.8.ht2l
steve.1.ht2l
...
steve.8.ht2l
和sereval RNAseq样本
flower_kevin_1.fastq.gz
flower_kevin_2.fastq.gz
flower_daniel_1.fastq.gz
flower_daniel_2.fastq.gz
我需要将所有rnaseq读数与每个基因组对齐。dariober建议更新:
workdir: "/path/to/aligned"
(HISAT2_INDEX_PREFIX,)=glob_wildcards("/path/to/index/{prefix}.1.ht2l")
(SAMPLES,)=glob_wildcards("/path/to/{sample}_1.fastq.gz")
print(HISAT2_INDEX_PREFIX)
print (SAMPLES)
rule all:
input:
expand("{prefix}.{sample}.bam", zip, prefix=HISAT2_INDEX_PREFIX, sample=SAMPLES)
rule hisat2:
input:
hisat2_index=expand("%s.{ix}.ht2l" % "/path/to/index/{prefix}", ix=range(1, 9), prefix = HISAT2_INDEX_PREFIX),
fastq1="/path/to/{sample}_1.fastq.gz",
fastq2="/path/to/{sample}_2.fastq.gz"
output:
bam = "{prefix}.{sample}.bam",
txt = "{prefix}.{sample}.txt",
log: "{prefix}.{sample}.snakemake_log.txt"
threads: 5
shell:
"/Tools/hisat2-2.1.0/hisat2 -p {threads} -x {/path/to/index/{wildcards.prefix}"
" -1 {input.fastq1} -2 {input.fastq2} --summary-file {output.txt} |"
"/Tools/samtools-1.9/samtools sort -@ {threads} -o {output.bam}"
我得到的问题是,在运行HISAT2时,一次将所有bob.1.ht2l:bob.8.ht2l和steve.1.ht2l:steve.8.ht2l都作为-x输入。虽然rna-seq应该分别定位在每个基因组上。错误在哪里?注意:我之前的问题:Snakemake: HISAT2 alignment of many RNAseq reads against many genomes
我认为您的困惑来自于事实,那就是hisat希望为索引文件提供prefix,而不是索引文件的所有列表。因此,请使用-x {input.hisat2_index}
代替-x /path/to/{wildcards.prefix}
(即索引文件列表)。
换句话说,输入hisat2_index=expand(...)
应该只在这些文件准备好之后才告诉snakemake启动此规则,但是您不能直接使用它们(当然,hisat确实会使用它们,但是您不会在命令行上传递它们。