how can I align several fastq files with hisat2 and get htseq-count

Question

[Learner]

Hello,

I have been strugelling to run Hisat2 code for alignment of several fastq files from human I have read as many post as I could online and I went through the protocol here https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5032908/

At the moment, i have the latest binary version of hisat2 in $HOME/bin When I do hisat2 in the terminal and shows that I don't import any input but it indicates that it is there and can be invoked

No index, query, or output file specified! HISAT2 version 2.1.0 by Daehwan Kim ([e-mail], www.ccb.jhu.edu/people/infphilo) Usage:

I have 10 fasta files on a folder named data on my desktop. I have put the three scripts attached to the protocl in that folder and the files are shown below

NIHMS816842-supplement-README.txt
RNAseq_shell_script.sh
SRR1_1.fastq.gz SRR9_1.fastq.gz SRR1_2.fastq.gz SRR9_2.fastq.gz SRR2_1.fastq.gz SRR10_1.fastq.gz SRR2_2.fastq.gz SRR10_2.fastq.gz SRR3_1.fastq.gz SRR3_2.fastq.gz SRR4_1.fastq.gz SRR4_2.fastq.gz SRR5_1.fastq.gz SRR5_2.fastq.gz SRR6_1.fastq.gz SRR6_2.fastq.gz SRR7_1.fastq.gz SRR7_2.fastq.gz rnaseq_ballgown.R SRR8_1.fastq.gz rnaseq_pipeline.config.sh SRR8_2.fastq.gz

Can you please let me know how I can run the code to align these fastq files with a Human gene reference ?

I did not change anything on any of the script which I believe that I must change the one in rnaseq_pipline.config.sh

At this moment, it is

NUMCPUS=8
HISAT2=$(which hisat2)
STRINGTIE=$(which stringtie)
SAMTOOLS=$(which samtools)
BASEDIR=$(pwd -P)/chrX_data
FASTQLOC="$BASEDIR/samples"
GENOMEIDX="$BASEDIR/indexes/chrX_tran"
GTFFILE="$BASEDIR/genes/chrX.gtf"
PHENODATA="$BASEDIR/geuvadis_phenodata.csv"
TEMPLOC="./tmp" #this will be relative to the output directory
reads1=(${FASTQLOC}/*_1.*)
reads1=("${reads1[@]##*/}")
reads2=("${reads1[@]/_1./_2.}")

All what I want is to get the htseq-count. txt files for each sample

based on this post Run mutltiple Fastq files on Hisat2 If I do the following

for f in `ls -1 *_1.fastq.gz | sed 's/_1.fastq.gz//' `
do
echo hisat2 -x mm10idx -1 ${f}_1.fastq.gz -2 ${f}_2.fastq.gz -S ${f}.bam
done

It will generate all the bam files but all empty. if I echo the do , it will be the output as an example

hisat2 -x mm10idx -1 SRR1_1.fastq.gz -2 SRR1_2.fastq.gz -S SRR1.bam
hisat2 -x mm10idx -1 SRR2_1.fastq.gz -2 SRR42_2.fastq.gz -S SRR2.bam
hisat2 -x mm10idx -1 SRR3_1.fastq.gz -2 SRR3_2.fastq.gz -S SRR3.bam
hisat2 -x mm10idx -1 SRR4_1.fastq.gz -2 SRR_2.fastq.gz -S SRR4.bam
hisat2 -x mm10idx -1 SRR5_1.fastq.gz -2 SRR5_2.fastq.gz -S SRR5.bam
hisat2 -x mm10idx -1 SRR6_1.fastq.gz -2 SRR6_2.fastq.gz -S SRR6.bam
hisat2 -x mm10idx -1 SRR7_1.fastq.gz -2 SRR7_2.fastq.gz -S SRR7.bam
hisat2 -x mm10idx -1 SRR8_1.fastq.gz -2 SRR8_2.fastq.gz -S SRR8.bam
hisat2 -x mm10idx -1 SRR9_1.fastq.gz -2 SRR9_2.fastq.gz -S SRR9.bam
hisat2 -x mm10idx -1 SRR10_1.fastq.gz -2 SRR10_2.fastq.gz -S SRR10.bam

On the other hand, where should I involve the reference genome? for example, I downloaded the hg19 from here https://ccb.jhu.edu/software/hisat2/index.shtml

[Hugo]

Maybe DEWE (Differential Workflow Executor for RNA-Seq) can be helpful to you: http://www.sing-group.org/dewe It provides a user friendly GUI to run HISAT2, convert aligments to bam and run htseq-count. Regards.

[Devon Ryan]

Better to use Galaxy, it's widely used and provides most of the common tools.

[Learner]

@Devon Ryan can you please tell me how to use galaxy for it? I know how to upload the data and how to run the hisat but some parameters I don't know. For example, single end or paired reads? I guess pair read? or input read as individual or paird? and the rest default values. but the question is how then get the htseq-count? because when I try to use htseq-count function, it needs gff file which I don't have

[Devon Ryan]

You have paired-end data, so you should select paired end. Run featureCounts on the resulting BAM files. It's in Galaxy too, as is DESeq2. In fact, there's training material available online (it's what we made to use when we train people here in Freiburg).

[Learner]

@Devon Ryan thank you . actually many people think of htseq-count rather than feature counts. I do't know which one is the best but I will try to follow your treating protocol and contact you if needed. Thank you again I believe Galaxy is amazing and easier to work with

[Learner]

@Devon Ryan can you please tell me how to get this gff for human ? can I use http://genome.ucsc.edu/cgi-bin/hgTables and then send it to galaxy or do you know a better way?

[Devon Ryan]

Never use GFF files or annotations from UCSC. Download a GTF file from Gencode or Ensembl.

[Learner]

@Devon Ryan can you please give me a link to the right gft? why would UCSC be a problem?

[Devon Ryan]

Here is the link to the Gencode, release 27 primary assembly GTF. UCSC annotations have historically both lagged behind those of Ensembl/Gencode and been of lower quality. They have also historically had a variety of odd issues, such as having multiple instances of the same gene on multiple strands/chromosomes all sharing the same ID...which makes no biological sense (they can have the same name, but they need different IDs since they're different loci). Ensembl/Gencode GTFs won't break tools, UCSC annotations have been known to do exactly that.

Answer 1

[Devon Ryan]

You're trying to do too many things at once, which is just going to cause confusion for someone new to the field. Start out with a single sample that manually type the command for. Once that works properly you can move on to multiple samples in a loop.

Since you apparently downloaded the hg19 prebuilt reference, you need to specify where it is with the -x option. In your examples, that was -x mm10idx, but presumably something like -x /some/path/to/the/downloaded/files/hg19 would make more sense.

[Learner]

@Devon Ryan I am trying to run this, I will accept your answer once I get it to run

[Learner]

@Devon Ryan that helped me thanks although I had to do a lot of work to get it to run :-D

[psschlogl]

@Devon Ryan Hi Devo I tried to run this bash script:

SAMPLES="x1..x9"

for SAMPLE in $SAMPLES; do
    hisat2 -p 4 --dta -x ${hisat2Ref} -U fastq/${SAMPLE}.fastq.gz | samtools view -bh - | samtools sort - bam/${SAMPLE}.bam | samtools index bam/${hisat2Ref}.bam
done

But I didn't get it working do you could give one idea why it didn't worked? I not experienced in bash scripting. Thanks

[Devon Ryan]

I strongly recommend that you not use bash scripting and instead use something more scalable and controllable like snakemake.

[psschlogl]

Thanks I will check it out.