Bulk RNA-Seq Analysis Series: Part 9 - HTSeq counts
Welcome to Part 9 of our Bulk RNA-Seq Analysis series! In this video, we'll be focusing on HTSeq for counting features, such as genes, in your aligned RNA-Seq data. HTSeq is a flexible and powerful tool that is commonly used in RNA-Seq analysis pipelines.
By the end of this tutorial, you will know how to:
🔹 Count features using HTSeq
🔹 Analyze the output for downstream applications
Script File Used for analysis
#### Run htseq-count for each BAM file
for bamfile in ./hisat2_bam_file/*_R_trimmomatic_sortmerna_sorted_bam; do
sample_name=$(basename $bamfile | cut -d'_' -f1)
htseq-count -f bam -r pos -s no -t exon -i gene_id $bamfile /home/shashi/rnaseq/Mus_musculus.GRCm39.112.gtf
./hisat2_count/gene_counts_${sample_name}.txt
done
#### Create header for combined counts file
header="GeneID"
for bamfile in ./hisat2_bam_file/*_R_trimmomatic_sortmerna_sorted_bam; do
sample_name=$(basename $bamfile | cut -d'_' -f1)
header="${header}\t${sample_name}"
done
echo -e $header
./hisat2_count/gene_counts_combined.txt
#### Extract gene IDs from the first count file
first_sample=$(ls ./hisat2_count/gene_counts_*.txt | head -n 1)
cut -f1 $first_sample
./hisat2_count/gene_ids.txt
Initialize the combined counts file with gene IDs
cp ./hisat2_count/gene_ids.txt ./hisat2_count/gene_counts_combined.txt
Combine counts from each file into the combined counts file
for countfile in ./hisat2_count/gene_counts_*.txt; do
cut -f2 $countfile
./hisat2_count/current_count.txt
paste ./hisat2_count/gene_counts_combined.txt ./hisat2_count/current_count.txt
./hisat2_count/gene_counts_tmp.txt
mv ./hisat2_count/gene_counts_tmp.txt ./hisat2_count/gene_counts_combined.txt
rm ./hisat2_count/current_count.txt
done
Clean up temporary file
rm ./hisat2_count/gene_ids.txt
Link for Mouse Genome GTF file
https://ftp.ensembl.org/pub/release-1...
Script description
for bamfile in ./hisat2_bam_file/*_R_trimmomatic_sortmerna_sorted_bam; do: Loops over all BAM files matching the pattern in the directory ./hisat2_bam_file/.
sample_name=$(basename $bamfile | cut -d'_' -f1): Extracts the base name of the BAM file (removes the directory path) and cuts the string at the first underscore (_) to get the sample name.
htseq-count -f bam -r pos -s no -t exon -i gene_id $bamfile /home/shashi/rnaseq/Mus_musculus.GRCm39.112.gtf
./hisat2_count/gene_counts_${sample_name}.txt: Runs htseq-count on each BAM file, specifying that the input format is BAM (-f bam), reads are sorted by position (-r pos), and the data is not strand-specific (-s no). The count results are written to a file named gene_counts_${sample_name}.txt in the ./hisat2_count/ directory.
header="GeneID": Initializes the header string with the text "GeneID", which will be the first column in the combined counts file.
for bamfile in ... ; do ... done: Loops over each BAM file again to add the sample names to the header.
sample_name=$(basename $bamfile | cut -d'_' -f1): Extracts the sample name as before.
header="${header}\t${sample_name}": Adds the sample name to the header, separated by a tab (\t).
echo -e $header
./hisat2_count/gene_counts_combined.txt: Writes the header to the combined counts file.
first_sample=$(ls ./hisat2_count/gene_counts_*.txt | head -n 1): Gets the first count file in the ./hisat2_count/ directory.
cut -f1 $first_sample
./hisat2_count/gene_ids.txt: Extracts the first column (Gene IDs) from the first count file and saves it to gene_ids.txt.
cp ./hisat2_count/gene_ids.txt ./hisat2_count/gene_counts_combined.txt: Copies the gene_ids.txt file to the combined counts file. This initializes the combined file with the Gene IDs.
for countfile in ./hisat2_count/gene_counts_*.txt; do ... done: Loops over all the individual count files.
cut -f2 $countfile
./hisat2_count/current_count.txt: Extracts the second column (the counts) from each count file and saves it as current_count.txt.
paste ./hisat2_count/gene_counts_combined.txt ./hisat2_count/current_count.txt
./hisat2_count/gene_counts_tmp.txt: Combines the existing combined file with the current counts column.
mv ./hisat2_count/gene_counts_tmp.txt ./hisat2_count/gene_counts_combined.txt: Moves the temporary combined file back to the main combined file.
rm ./hisat2_count/current_count.txt: Deletes the temporary count file.
rm ./hisat2_count/gene_ids.txt: Removes the temporary gene_ids.txt file as it's no longer needed.
This script takes multiple BAM files, counts the reads per gene using htseq-count, and then combines these counts into a single file where each column corresponds to a sample, and each row corresponds to a gene. The paste command is used to concatenate the counts for each sample into a final combined table.
#Learn Innovatively with Me