Compare commits
No commits in common. "master" and "step3/star-alignment-snippet" have entirely different histories.
master
...
step3/star
|
|
@ -1,12 +1,5 @@
|
|||
#!/bin/bash
|
||||
R1="sample1_R1.fastq.gz"
|
||||
R2="sample1_R2.fastq.gz"
|
||||
# The command to do a FastQC on a fastq file is
|
||||
file="file_to_analyse.fq.gz"
|
||||
fastqc_out="./path/to/fastqc/output/dir"
|
||||
|
||||
PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
|
||||
FASTQC_OUT="${PROJECT_DIRECTORY}/step1/"
|
||||
mkdir -p "${FASTQC_OUT}"
|
||||
|
||||
# Run FastQC on paired-end data.
|
||||
fastqc \
|
||||
-o "${FASTQC_OUT}" \
|
||||
"${R1}" "${R2}"
|
||||
fastqc -o "$fastqc_out" "$file"
|
||||
|
|
|
|||
|
|
@ -1,12 +1,9 @@
|
|||
#!/bin/bash
|
||||
#
|
||||
# Reference: http://www.usadellab.org/cms/?page=trimmomatic
|
||||
|
||||
# reference: http://www.usadellab.org/cms/?page=trimmomatic
|
||||
|
||||
|
||||
module load Trimmomatic
|
||||
PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
|
||||
FASTQ_OUT="${PROJECT_DIRECTORY}/step2/"
|
||||
mkdir -p "${FASTQ_OUT}"
|
||||
|
||||
|
||||
# Adapters can be found at
|
||||
|
|
@ -14,26 +11,16 @@ mkdir -p "${FASTQ_OUT}"
|
|||
# But should be verified with FastQC, or in another way.
|
||||
|
||||
|
||||
# Trimmomatic example Paired end data.
|
||||
#
|
||||
# Flags:
|
||||
# - ILLUMINACLIP: Cut adapter and other illumina-specific sequences from the
|
||||
# read.
|
||||
# - SLIDINGWINDOW: Perform a sliding window trimming, cutting once the average
|
||||
# quality within the window falls below a threshold.
|
||||
# - LEADING: Cut bases off the start of a read, if below a threshold quality.
|
||||
# - TRAILING: Cut bases off the end of a read, if below a threshold quality.
|
||||
# - HEADCROP: Cut the specified number of bases from the start of the read.
|
||||
# - MINLEN: Drop the read if it is below a specified length.
|
||||
# (Example) Paired end
|
||||
java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar PE \
|
||||
-phred33 \
|
||||
sample1_R1.fastq.gz \
|
||||
sample1_R2.fastq.gz \
|
||||
"${FASTQ_OUT}/sample1_R1_paired.fastq.gz" \
|
||||
"${FASTQ_OUT}/sample1_R1_unpaired.fastq.gz" \
|
||||
"${FASTQ_OUT}/sample1_R2_paired.fastq.gz" \
|
||||
"${FASTQ_OUT}/sample1_R2_unpaired.fastq.gz" \
|
||||
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
|
||||
input_forward.fq.gz \
|
||||
input_reverse.fq.gz \
|
||||
output_forward_paired.fq.gz \
|
||||
output_forward_unpaired.fq.gz \
|
||||
output_reverse_paired.fq.gz \
|
||||
output_reverse_unpaired.fq.gz \
|
||||
ILLUMINACLIP: TruSeq3-PE.fa:2:30:10 \
|
||||
LEADING:3 \
|
||||
TRAILING:3 \
|
||||
SLIDINGWINDOW:4:25 \
|
||||
|
|
@ -41,11 +28,11 @@ java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar PE \
|
|||
MINLEN:50
|
||||
|
||||
|
||||
# Example single end data.
|
||||
# (Example) Single end
|
||||
java -jar $EBROOTTRIMMOMATIC/trimmomatic.jar SE \
|
||||
-phred33 \
|
||||
sample1.fastq.gz \
|
||||
output.fastq.gz \
|
||||
input.fq.gz \
|
||||
output.fq.gz \
|
||||
ILLUMINACLIP:TruSeq3-SE:2:30:10 \
|
||||
LEADING:3 \
|
||||
TRAILING:3 \
|
||||
|
|
|
|||
|
|
@ -2,42 +2,32 @@
|
|||
#
|
||||
# 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"
|
||||
STORAGE="/groups/umcg-griac/tmp04/rawdata/$(whoami)/step3"
|
||||
# Store genome index in this location:.
|
||||
GENOME_INDEX="${STORAGE}/genome_index"
|
||||
mkdir -p "${GENOME_INDEX}"
|
||||
# Store the generated `Aligned.sortedByCoord.out.bam` in this dir.
|
||||
ALIGNMENT_OUTPUT="${STORAGE}/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.
|
||||
# 1) Generate genome index.
|
||||
#
|
||||
# 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
|
||||
# - We're assuming a read size of 100 bp (--sjdbOverhang 100). 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.
|
||||
# - We're using gzip compressed reference data (--readFilesCommand zcat), i.e.,
|
||||
# .gtf.gz and fa.gz. If not, you can remove the `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"
|
||||
# Storage location reference data (in this case on calculon).
|
||||
REFERENCE_DATA="/groups/umcg-griac/prm02/rawdata/reference/genome"
|
||||
GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf.gz"
|
||||
FASTA_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"
|
||||
|
||||
STAR \
|
||||
--runThreadN 8 \
|
||||
--runMode genomeGenerate \
|
||||
--readFilesCommand zcat \
|
||||
--sjdbOverhang 100 \
|
||||
--genomeFastaFiles ${FASTA_FILE} \
|
||||
--sjdbGTFfile ${GTF_FILE} \
|
||||
|
|
@ -50,9 +40,9 @@ STAR \
|
|||
# - 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"
|
||||
# THe compressed paired-end FastQ's that we are aligning.
|
||||
R1="sample1_R1.fastq.gz"
|
||||
R2="sample1_R2.fastq.gz"
|
||||
|
||||
STAR \
|
||||
--runThreadN 8 \
|
||||
|
|
@ -60,4 +50,4 @@ STAR \
|
|||
--readFilesIn "${R1}" "${R2}" \
|
||||
--outSAMtype BAM SortedByCoordinate \
|
||||
--genomeDir ${GENOME_INDEX} \
|
||||
--outFileNamePrefix "${ALIGNMENT_OUTPUT}/sample1_"
|
||||
--outFileNamePrefix "${ALIGNMENT_OUTPUT}"
|
||||
|
|
|
|||
|
|
@ -1,5 +0,0 @@
|
|||
#!/bin/bash
|
||||
module load multiqc
|
||||
|
||||
PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
|
||||
multiqc "${PROJECT_DIRECTORY}"
|
||||
|
|
@ -1,28 +0,0 @@
|
|||
#!/bin/bash
|
||||
PROJECT_DIRECTORY="/groups/umcg-griac/tmp01/rawdata/$(whoami)/rnaseq"
|
||||
COUNT_OUTPUT="${PROJECT_DIRECTORY}/step5"
|
||||
mkdir -p "${COUNT_OUTPUT}"
|
||||
|
||||
# Storage location of annotation on Gearshift.
|
||||
REFERENCE_DATA="/groups/umcg-griac/prm03/rawdata/reference/genome"
|
||||
GTF_FILE="${REFERENCE_DATA}/Homo_sapiens.GRCh38.100.gtf"
|
||||
|
||||
# Where our alignment file was stored.
|
||||
BAM="${PROJECT_DIRECTORY}/step3/alignment/sample1_Aligned.sortedByCoord.out.bam"
|
||||
|
||||
# Compute counts using htseq-count.
|
||||
#
|
||||
# N.B.:
|
||||
# - If you are processing multiple files, consider using the `--nprocesses` flag
|
||||
# to distribute computation of the files to different cores.
|
||||
# - The BAM file must be position sorted. If you used STAR with the
|
||||
# `SortedByCoordinate` option you should be okay. If not, sort your BAM file
|
||||
# using `samtools sort`.
|
||||
# - By default, strand aware library preparation is assumed. If not, specify the
|
||||
# `--stranded` flag.
|
||||
htseq-count \
|
||||
--order pos \
|
||||
-f bam \
|
||||
${BAM} \
|
||||
${GTF_FILE} \
|
||||
> ${COUNT_OUTPUT}/counts.txt
|
||||
|
|
@ -18,7 +18,7 @@ do.scale = FALSE
|
|||
|
||||
|
||||
# The analysis
|
||||
# We use prcomp to calculate the PCAs. Afterwards you should plot the results.
|
||||
# We ise prcomp to calculate the PCAs. Afterwards you should plot the results.
|
||||
norm.expr.data <- expression.data %>%
|
||||
tibble::column_to_rownames("Gene")
|
||||
norm.expr.data <- norm.expr.data[rowSums(norm.expr.data) >= 10,] %>%
|
||||
|
|
@ -37,8 +37,8 @@ norm.expr.data.pcs <- norm.expr.data %>%
|
|||
)
|
||||
|
||||
# Write summary of PCAs to files
|
||||
pcs.summary <- summary(norm.expr.data.pcs)
|
||||
pcs.summary$importance %>%
|
||||
pcs.summery <- summary(norm.expr.data.pcs)
|
||||
pcs.summery$importance %>%
|
||||
t() %>%
|
||||
as.data.frame() %>%
|
||||
tibble::rownames_to_column("PC.name") %>%
|
||||
|
|
@ -46,7 +46,7 @@ pcs.summary$importance %>%
|
|||
file.path(results.dir.pca, "importance.csv")
|
||||
)
|
||||
|
||||
pcs.summary$x %>%
|
||||
pcs.summery$x %>%
|
||||
t() %>%
|
||||
as.data.frame() %>%
|
||||
tibble::rownames_to_column("ensembl.id") %>%
|
||||
|
|
@ -54,7 +54,7 @@ pcs.summary$x %>%
|
|||
file.path(results.dir.pca, "values.csv")
|
||||
)
|
||||
|
||||
pcs.summary$rotation %>%
|
||||
pcs.summery$rotation %>%
|
||||
t() %>%
|
||||
as.data.frame() %>%
|
||||
tibble::rownames_to_column("sample.id") %>%
|
||||
|
|
@ -63,15 +63,15 @@ pcs.summary$rotation %>%
|
|||
)
|
||||
|
||||
data.frame(
|
||||
rownames = names(pcs.summary$center),
|
||||
center = pcs.summary$center,
|
||||
scale = pcs.summary$scale
|
||||
rownames = names(pcs.summery$center),
|
||||
center = pcs.summery$center,
|
||||
scale = pcs.summery$scale
|
||||
) %>%
|
||||
readr::write_csv(
|
||||
file.path(results.dir.pca, "rest.csv")
|
||||
)
|
||||
|
||||
# Not saved: pcs.summary$sdev,
|
||||
# Not saved: pcs.summery$sdev,
|
||||
|
||||
|
||||
# Next thing to do:
|
||||
|
|
|
|||
|
|
@ -1,61 +0,0 @@
|
|||
#!/usr/bin/perl -w
|
||||
use strict;
|
||||
use Parallel::ForkManager;
|
||||
# this script creats one file with UMI unique reads and one with UMI duplicated reads
|
||||
# as input you need aligned sorted by coordinate bam file
|
||||
my @torun = ();
|
||||
foreach my $file ( <*Aligned.sortedByCoord.out.bam> ) {
|
||||
push @torun, $file;
|
||||
}
|
||||
|
||||
my $pm = Parallel::ForkManager->new( 12 );
|
||||
foreach my $file ( @torun ) {
|
||||
my $sample = $file;
|
||||
$sample =~ s/Aligned\.sortedByCoord\.out\.bam$//;
|
||||
next if -s $sample.'_uniq.bam';
|
||||
warn "Parsing $sample\n";
|
||||
$pm->start and next;
|
||||
my %seen = ();
|
||||
my %duplicates = ();
|
||||
my ( $uniqs, $dups ) = ( 0, 0 );
|
||||
open F, 'samtools view -h '.$file.' |';
|
||||
open F1, '| samtools view -bS - >'.$sample.'_uniq.bam';
|
||||
open F2, '| samtools view -bS - >'.$sample.'_dups.bam';
|
||||
while ( <F> ) {
|
||||
if ( m/^\@/ ) {
|
||||
print F1;
|
||||
print F2;
|
||||
next;
|
||||
}
|
||||
my ( $id, $flag, $chr, $pos, $mapq, $cigar, $chr2, $pos2, $tlen ) = split /\t/;
|
||||
next if $flag & 256 or $flag & 512 or $flag & 1024; #skip if the read is not primary alignment/read fails platform/vendor quality checks/read is PCR or optical duplicate
|
||||
# foreach ( 256, 512, 1024 ) { $flag-=$_ if $flag&$_ }
|
||||
my ( $bc ) = $id =~ m/\:([GATCN\d]+)$/; #extract UMI barcode
|
||||
my $uniq = join( ':', $chr, $pos, $flag, $tlen, $bc );
|
||||
my $pos_ = $pos-1;
|
||||
while ( $cigar =~ m/(\d+)([SHMDIN=])/g ) {
|
||||
$pos_+=$1 if $2 eq 'M' or $2 eq '=' or $2 eq 'D' or $2 eq 'N'; # find position for minus strand
|
||||
}
|
||||
my $uniq2 = join( ':', $chr, $pos_, $flag, $tlen, $bc );
|
||||
if ( exists($duplicates{$id}) or # already marked as duplicate
|
||||
( not($flag&16) and ++$seen{ $uniq } > 1 ) or # plus strand
|
||||
( $flag&16 and ++$seen{ $uniq2 } > 1 ) #minus strand
|
||||
) {
|
||||
print F2;
|
||||
$dups++;
|
||||
$duplicates{$id} = 1;
|
||||
}
|
||||
else {
|
||||
print F1;
|
||||
$uniqs++;
|
||||
}
|
||||
}
|
||||
close F;
|
||||
close F1;
|
||||
close F2;
|
||||
warn "Unique: $uniqs (", 100*$uniqs/($uniqs+$dups), ")\n";
|
||||
system( 'samtools', 'index', $sample.'_uniq.bam' );
|
||||
# last;
|
||||
$pm->finish;
|
||||
}
|
||||
$pm->wait_all_children;
|
||||
Loading…
Reference in New Issue