-
Notifications
You must be signed in to change notification settings - Fork 1
3. Bulk RNAseq
Karla Lozano Gonzalez edited this page Jun 16, 2023
·
4 revisions
cd /path/to/genome_dir
mkdir -p Hisat2_index
# The name of the index files will be as genome.x.ht# and will be saved in the Hisat2_index directory
hisat2-build -p 38 genome.fa Hisat2_index/genome
# Generate file with known splice junctions to be fed to HISAT2
hisat2_extract_exons.py genome.gtf > genome_known_splice_sites.txtcd /PATH/TO/PROJECT/outputs
mkdir -p Hisat2_SAMs
# parse fastq data into arrays
for fwd in PATH/TO/PROJECT/originals/FASTQ/*.R1.fastq.gz;do
rev=${fwd/R1/R2} # Change the pattern and replacement accordingly
name=$( basename $fwd )
name=${name%.R1.fastq.gz} # Change pattern accordingly
hisat2 -p 20 --dta --known-splicesite-infile /path/to/genome_known_splice_sites.txt -x /path/to/Hisat2_index/genome \
-1 $fwd -2 $rev -S Hisat2_SAMs/${name}.sam --no-discordant
doneFor Single-end reads, the same code can be used, but replace -1 with -U and exclude the pair read.
Because HISAT outputs SAM files, we will have to convert to a more portable BAM format.
# sort and convert .sam files into BAM format
cd /PATH/TO/PROJECTS/output
mkdir Hisat2_sorted_BAMs
for sam in Hisat2_SAMs/*.sam; do
name=$( basename $sam )
name=${name%.sam}
samtools view -@ 20 -Su $sam | samtools sort -@ 20 -o Hisat2_sorted_BAMs/$name.bam
done
# index bam files
cd Hisat2_sorted_BAMs
for filename in *.bam; do
samtools index -@ 40 $filename $filename.bai
doneWe use StringTie2 to construct our own custom transcriptome containing de novo transcripts.
cd /PATH/TO/PROJECTS/output
mkdir -p Stringtie_gtf
for bam in Hisat2_sorted_BAMs/*.bam; do
name=$( basename $bam )
name=${name%.bam}
stringtie $bam -p 20 -o Stringtie_gtf/$name.gtf -G /path/to/referenceGTF
done
Common lab SOPs:
Bioinformatics-related:
Image analyses-related:
Programming-related: