Hello,
The general command for Hisat2 alignment is:
hisat2 --threads 8 -t -p 8 --summary-file summary_txt -x /hisat2_index/index -1 read1.fastq.gz -2 read2.fastq.gz -S aligned.sam
I use the following script to run my Hisat2 pipelines. The bam files I am getting are fine. My question is how can I get the summary file for alignment too? How to frame the output to capture the summary.txt files for each alignment?
HISAT2_INDEX_PREFIX = "/hisat2_index/index"
SAMPLES, *_=glob_wildcards('/samples/{sample}.read1.fastq.gz')
rule align_all_samples:
input: expand("samples/{sample}.bam", sample=SAMPLES)
rule align_hisat2:
input:
hisat2_index=expand(f"{HISAT2_INDEX_PREFIX}.{{ix}}.ht2", ix=range(1, 9)),
fastq1="/samples/{sample}.read1.fastq.gz",
fastq2="/samples/{sample}.read2.fastq.gz",
output:
bam = "/samples/align_hisat2_update1/{sample}.bam",
log: "/samples/snakemake_log.txt"
threads: 8
shell:
"hisat2 -p {threads} -x {HISAT2_INDEX_PREFIX}"
" -1 {input.fastq1} -2 {input.fastq2} |"
"samtools sort -@ {threads} -o {output.bam}"
I tried taking another output variable like summary file under the outputs section, but that gives me error! Can I club two different types of output files in the rule_all?
Thanks, Payal
•
written
by
Payal
The stderr method suggested by Devon should work, but it is cleaner to just use the built-in --summary-file option. Just modify the "shell" part of your script into something like this:
"hisat2 -p {threads} --summary-file {sample}.summary.txt -x {HISAT2_INDEX_PREFIX} -1 {input.fastq1} -2 {input.fastq2} | samtools sort -@ {threads} -o {output.bam}"
Thank you Carlo. I tried this but it gives me
NameError:
The name 'sample' is unknown in this context. Please make sure that you defined that variable. Also note that braces not used for variable access have to be escaped by repeating them, i.e. {{print $1}}
How should I fix this?
I'm not really familiar with snakemake, so I'm not sure, but you can perhaps try this:
rule align_hisat2:
input:
hisat2_index=expand(f"{HISAT2_INDEX_PREFIX}.{{ix}}.ht2", ix=range(1, 9)),
fastq1="/samples/{sample}.read1.fastq.gz",
fastq2="/samples/{sample}.read2.fastq.gz",
output:
bam = "/samples/align_hisat2_update1/{sample}.bam",
txt = "/samples/align_hisat2_update1/{sample}.txt",
log: "/samples/snakemake_log.txt"
threads: 8
shell:
"hisat2 -p {threads} -x {HISAT2_INDEX_PREFIX}"
" -1 {input.fastq1} -2 {input.fastq2} --summary-file {output.txt} |"
"samtools sort -@ {threads} -o {output.bam}"
The alignment summary should (hopefully) be saved in the "/samples/align_hisat2_update1/{sample}.txt" file.
Thank you very much Carlo. I am also new to snakemake, so I test and periodically try to improve my script.
Actually, I just tested my scirpt using the same logic that you have described here. It works perfectly now. I was going to post my updated script, coincidentally you have already posted it with the same logic.
Cheers! Payal
Hi, may I know how do you interpret the output using HISAT2? I use the data from the public paper (Pertea et al. - 2016 - Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown) to learn RNAseq analysis.
The command is 'hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR188044_chrX_1.fastq.gz -2 chrX_data/samples/ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam'. and the log is shown below.
1321477 reads; of these: 1321477 (100.00%) were paired; of these:
121980 (9.23%) aligned concordantly 0 times 1183228 (89.54%) aligned concordantly exactly 1 time 16269 (1.23%) aligned concordantly >1 times ---- 121980 pairs aligned concordantly 0 times; of these: 4166 (3.42%) aligned discordantly 1 time ---- 117814 pairs aligned 0 times concordantly or discordantly; of these: 235628 mates make up the pairs; of these: 119930 (50.90%) aligned 0 times 113038 (47.97%) aligned exactly 1 time 2660 (1.13%) aligned >1 times95.46% overall alignment rate
Sample id Uniquely aligned reads (%) Multimapped reads (%) Overall alignment rate (%)
ERR188044 73.1 21.9 95.0
Please post this as a new question and mention which part you do not understand.
Log in to answer this question.
Isn't that printed to
stderr? So just redirect that to a new file.I am sorry, didn't understand. Can you please elaborate?
Please google what
stderris in the Unix world and how you can save things printed to screen to a file. https://askubuntu.com/questions/625224/how-to-redirect-stderr-to-a-fileThank you Devon and AT.
I see what you are saying. I will try this. Will this output all sample summaries in one summary file like a log file or individual samples will output their individual summary report?