64 lines
2.2 KiB
Bash
64 lines
2.2 KiB
Bash
#!/bin/bash
|
|
#
|
|
# Align reads against reference genome.
|
|
|
|
REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
|
|
|
|
PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
|
|
# Store the generated `sample1_Aligned.sortedByCoord.out.bam` in this dir.
|
|
ALIGNMENT_OUTPUT="${PROJECT_DIRECTORY}/step3/alignment"
|
|
mkdir -p "${ALIGNMENT_OUTPUT}"
|
|
|
|
# 1) Generate genome index (optional).
|
|
#
|
|
# Depending on your read size, reference genome and annotation, you may need to
|
|
# generate a new genome index. In most cases, this is not necessary and you can
|
|
# directly use the pre-build genome index from the cluster:
|
|
#
|
|
# GENOME_INDEX="${REFERENCE_DATA}/index_GRCh38_gtf100_overhang100"
|
|
#
|
|
# and ignore the first STAR command below.
|
|
#
|
|
# N.B.:
|
|
# - We're assuming a read size of 100 bp (--sjdbOverhang 100). An alternative
|
|
# cut-off is 150, for low-input methods. In general, refer back to the
|
|
# previous quality control steps if you are unsure about the size. In case of
|
|
# reads of varying length, the ideal value is max(ReadLength)-1.
|
|
# - If you're using gzip compressed reference data, i.e., .gtf.gz and fa.gz,
|
|
# pass the `--readFilesCommand zcat` flag.
|
|
|
|
# Store created genome index in this location:.
|
|
GENOME_INDEX="${PROJECT_DIRECTORY}/step3/genome_index"
|
|
mkdir -p "${GENOME_INDEX}"
|
|
|
|
# Storage location reference data on Gearshift.
|
|
GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf"
|
|
FASTA_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.dna.primary_assembly.fa"
|
|
|
|
STAR \
|
|
--runThreadN 8 \
|
|
--runMode genomeGenerate \
|
|
--sjdbOverhang 100 \
|
|
--genomeFastaFiles ${FASTA_FILE} \
|
|
--sjdbGTFfile ${GTF_FILE} \
|
|
--genomeDir ${GENOME_INDEX}
|
|
|
|
|
|
# 2) Do the actual alignment.
|
|
#
|
|
# N.B.:
|
|
# - We are assuming paired-end, gzip compressed (--readFilesCommand zcat) FastQ
|
|
# files.
|
|
|
|
# The compressed, paired-end, FastQ's after trimming (step 2).
|
|
R1="${PROJECT_DIRECTORY}/step2/sample1_R1_paired.fastq.gz"
|
|
R2="${PROJECT_DIRECTORY}/step2/sample1_R2_paired.fastq.gz"
|
|
|
|
STAR \
|
|
--runThreadN 8 \
|
|
--readFilesCommand zcat \
|
|
--readFilesIn "${R1}" "${R2}" \
|
|
--outSAMtype BAM SortedByCoordinate \
|
|
--genomeDir ${GENOME_INDEX} \
|
|
--outFileNamePrefix "${ALIGNMENT_OUTPUT}/sample1_"
|