-
Notifications
You must be signed in to change notification settings - Fork 0
Code
Mathilda Stigenberg edited this page May 15, 2022
·
35 revisions
In order to run the program FastQC, the modules bioinfo-tools and FastQC/0.11.9 were loaded, as can be seen below.
module load bioinfo-tools FastQC/0.11.9
Then, when standing in the appropriate directory for these quality analyses, these commands below were run in order to obtain quality control for the raw data.
fastqc -o /home/mathi/genome_analysis/durianproject/analyses/1_fastqc_illumina/ -t 2 *.fastq.gz
fastqc -o /home/mathi/genome_analysis/durianproject/analyses/2_fastqc_pacbio/ -t 2 *.fq.gz
fastqc -o /home/mathi/genome_analysis/durianproject/analyses/3_fastqc_transcriptome/trimmed/ -t 2 *.fastq.gz
fastqc -o /home/mathi/genome_analysis/durianproject/analyses/3_fastqc_transcriptome/untrimmed/ -t 2 *.fastq.gz
Following the trimming of the untrimmed transcriptome with trimmomatic, a quality control of the recently trimmed trancriptome was done with the command below.
fastqc -o /home/mathi/genome_analysis/durianproject/analyses/4_fastqc/4_fastqc_transcriptome_trimmed/ -t 2 *.fq.gz
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 4
#SBATCH -t 19:00:00
#SBATCH -J 2_assembly_pacbio
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load canu/2.2
# Your commands
canu -d /home/mathi/genome_analysis/durianproject/analyses/2_assembly -p ass_pb genomeSize=24m useGrid=False corThreads=4 -pacbio /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/pacbio_data/SRR6037732_scaffold_10.fq.gz
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 2
#SBATCH -t 2:00:00
#SBATCH -J 2_trim_transcript
#SBATCH -o trimming_transcriptome.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load trimmomatic/0.39
# Your commands
trimmomatic PE -threads 2 /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/untrimmed/SRR6040095_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/untrimmed/SRR6040095_scaffold_10.2.fastq.gz -baseout /home/mathi/genome_analysis/durianproject/analyses/3_trim/3_trim_transcript/SRR6040095_scaffold_10.fq.gz ILLUMINACLIP:/sw/bioinfo/trimmomatic/0.39/snowy/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 2
#SBATCH -t 1:30:00
#SBATCH -J 5_mapping_bwa
#SBATCH -o mapping_bwa.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load bwa/0.7.17
module load samtools/1.14
# Your commands
bwa index /home/mathi/genome_analysis/durianproject/analyses/2_assembly/ass_pb.contigs.fasta
bwa mem -t 2 /home/mathi/genome_analysis/durianproject/analyses/2_assembly/ass_pb.contigs.fasta /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/illumina_data/SRR6058604_scaffold_10.1P.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/illumina_data/SRR6058604_scaffold_10.2P.fastq.gz > /home/mathi/genome_analysis/durianproject/analyses/5_mapping_bwa/aln_seq_bwa.sam
samtools sort /home/mathi/genome_analysis/durianproject/analyses/5_mapping_bwa/aln_seq_bwa.sam -o /home/mathi/genome_analysis/durianproject/analyses/5_mapping_bwa/aln_seq_bwa.bam
samtools index /home/mathi/genome_analysis/durianproject/analyses/5_mapping_bwa/aln_seq_bwa.bam /home/mathi/genome_analysis/durianproject/analyses/5_mapping_bwa/aln_seq_bwa.bai
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 2
#SBATCH -t 3:00:00
#SBATCH -J 6_pilon
#SBATCH -o pilon.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load Pilon/1.24
# Your commands
java -jar $PILON_HOME/pilon.jar --genome /home/mathi/genome_analysis/durianproject/analyses/2_assembly/ass_pb.contigs.fasta --bam /home/mathi/genome_analysis/durianproject/analyses/5_mapping_bwa/aln_seq_bwa.bam --outdir /home/mathi/genome_analysis/durianproject/analyses/6_improvement_pilon --threads 2
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 2
#SBATCH -t 15:00
#SBATCH -J 7_quast
#SBATCH -o quast.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load quast/5.0.2
# Your commands
python /sw/bioinfo/quast/5.0.2/snowy/bin/quast.py -t 2 -o /home/mathi/genome_analysis/durianproject/analyses/7_quality_assessment_quast -e /home/mathi/genome_analysis/durianproject/analyses/6_improvement_pilon/pilon.fasta
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 2
#SBATCH -t 05:00
#SBATCH -J 8_multiqc
#SBATCH -o multiqc.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load MultiQC/1.12
# Your commands
multiqc /home/mathi/genome_analysis/durianproject/analyses/1_fastqc/1_fastqc_illumina/*10* -o /home/mathi/genome_analysis/durianproject/analyses/8_multiqc/8_multiqc_illumina
multiqc /home/mathi/genome_analysis/durianproject/analyses/1_fastqc/1_fastqc_transcriptome/trimmed/*10* -o /home/mathi/genome_analysis/durianproject/analyses/8_multiqc/8_multiqc_transcriptome/8_trimmed
multiqc /home/mathi/genome_analysis/durianproject/analyses/1_fastqc/1_fastqc_transcriptome/untrimmed/*10* -o /home/mathi/genome_analysis/durianproject/analyses/8_multiqc/8_multiqc_transcriptome/8_untrimmed
multiqc /home/mathi/genome_analysis/durianproject/analyses/4_fastqc/4_fastqc_transcriptome_trimmed/*10* -o /home/mathi/genome_analysis/durianproject/analyses/8_multiqc/8_multiqc_transcriptome/8_untrimmed_is_trimmed
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 2
#SBATCH -t 1:00:00
#SBATCH -J 9_repeatmasker
#SBATCH -o repeatmasker.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load RepeatMasker/4.1.2-p1
# Your commands
RepeatMasker -pa 2 -species malvaceae -dir /home/mathi/genome_analysis/durianproject/analyses/9_repeatmasker /home/mathi/genome_analysis/durianproject/analyses/6_improvement_pilon/pilon.fasta
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 8
#SBATCH -t 30:00
#SBATCH -J 10_star
#SBATCH -o star.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load star/2.7.9a
module load samtools/1.14
# Your commands
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --genomeFastaFiles /home/mathi/genome_analysis/durianproject/analyses/9_repeatmasker/pilon.fasta.masked
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040092_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040092_scaffold_10.2.fastq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040092_scaffold_10_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040093_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040093_scaffold_10.2.fastq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040093_scaffold_10_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040094_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040094_scaffold_10.2.fastq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040094_scaffold_10_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040096_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040096_scaffold_10.2.fastq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040096_scaffold_10_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040097_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6040097_scaffold_10.2.fastq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040097_scaffold_10_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6156066_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6156066_scaffold_10.2.fastq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6156066_scaffold_10_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6156067_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6156067_scaffold_10.2.fastq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6156067_scaffold_10_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6156069_scaffold_10.1.fastq.gz /home/mathi/genome_analysis/durianproject/4_Tean_Teh_2017/transcriptome/trimmed/SRR6156069_scaffold_10.2.fastq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6156069_scaffold_10_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/analyses/3_trim/3_trim_transcript/SRR6040095_scaffold_10_1P.fq.gz /home/mathi/genome_analysis/durianproject/analyses/3_trim/3_trim_transcript/SRR6040095_scaffold_10_2P.fq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040095_scaffold_10_P_mapped --outSAMtype BAM SortedByCoordinate
STAR --runThreadN 8 --genomeDir /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/10_indexed_genome --readFilesIn /home/mathi/genome_analysis/durianproject/analyses/3_trim/3_trim_transcript/SRR6040095_scaffold_10_1U.fq.gz /home/mathi/genome_analysis/durianproject/analyses/3_trim/3_trim_transcript/SRR6040095_scaffold_10_2U.fq.gz --readFilesCommand zcat --outFileNamePrefix /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040095_scaffold_10_U_mapped --outSAMtype BAM SortedByCoordinate
Then, index files for the bam files were made with samtools index in the interactive session. The command can be seen below.
module load bioinfo-tools
module load samtools
for x in /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/*.bam
do
samtools index $x &
done
module load bioinfo-tools
module load augustus/20200915-ed0cb90
In order to copy the augustus configuration file to the desired directory (/home/mathi/genome_analysis/durianproject/analyses/11_braker_annotation), the modules above were loaded and then the command below was run.
source $AUGUSTUS_CONFIG_COPY
Then, the following commands below were run in order to have writing access to the configuration file and to copy the GeneMark key to the home directory.
chmod a+w -R /home/mathi/genome_analysis/durianproject/analyses/11_braker_annotation/augustus_config/species/
cp -vf /sw/bioinfo/GeneMark/4.33-es/snowy/gm_key $HOME/.gm_key
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 8
#SBATCH -t 2:00:00
#SBATCH -J 11_braker
#SBATCH -o braker.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load samtools/1.8
module load braker/2.1.1_Perl5.24.1
module load augustus/3.2.3_Perl5.24.1
module load bamtools/2.5.1
module load blast/2.9.0+
module load GenomeThreader/1.7.0
module load GeneMark/4.33-es_Perl5.24.1
# Your commands
cd /home/mathi/genome_analysis/durianproject/analyses/11_braker_annotation
/sw/bioinfo/braker/2.1.1/snowy/braker.pl --species durio_zibethinus --genome /home/mathi/genome_analysis/durianproject/analyses/9_repeatmasker/pilon.fasta.masked --bam /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040092_scaffold_10_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040093_scaffold_10_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040094_scaffold_10_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040095_scaffold_10_P_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040095_scaffold_10_U_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040096_scaffold_10_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6040097_scaffold_10_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6156066_scaffold_10_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6156067_scaffold_10_mappedAligned.sortedByCoord.out.bam,/home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/SRR6156069_scaffold_10_mappedAligned.sortedByCoord.out.bam --cores 8 --useexisting --AUGUSTUS_CONFIG_PATH=/home/mathi/genome_analysis/durianproject/analyses/11_braker_annotation/augustus_config --AUGUSTUS_BIN_PATH=/sw/bioinfo/augustus/3.4.0/snowy/bin --AUGUSTUS_SCRIPTS_PATH=/sw/bioinfo/augustus/3.4.0/snowy/scripts --GENEMARK_PATH=/sw/bioinfo/GeneMark/4.33-es/snowy --gff3
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 4
#SBATCH -t 12:00:00
#SBATCH -J 12_eggnogmapper
#SBATCH -o eggnogmapper.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load eggNOG-mapper/2.1.5
# Your commands
/sw/bioinfo/eggNOG-mapper/2.1.5/snowy/bin/emapper.py -i /home/mathi/genome_analysis/durianproject/analyses/9_repeatmasker/pilon.fasta.masked --itype genome -o functional_annotation --output_dir /home/mathi/genome_analysis/durianproject/analyses/12_eggnog_annotation --cpu 4
#!/bin/bash -l
#SBATCH -A uppmax2022-2-5
#SBATCH -M snowy
#SBATCH -p core
#SBATCH -n 8
#SBATCH -t 5:00:00
#SBATCH -J 13_htseq
#SBATCH -o htseq.output
#SBATCH --mail-type=ALL
#SBATCH --mail-user mathilda.stigenberg.5156@student.uu.se
# Load modules
module load bioinfo-tools
module load htseq/0.12.4
# Your commands
cd /home/mathi/genome_analysis/durianproject/analyses/13_htseq_countfeatures
for x in /home/mathi/genome_analysis/durianproject/analyses/10_star_mapping/*.bam
do
htseq-count -f bam -r pos -i ID -n 8 -o count -c ${x%%.*}_counted $x /home/mathi/genome_analysis/durianproject/analyses/11_braker_annotation/braker/durio_zibethinus/augustus.hints.gff3
done
# load packages
library("DESeq2")
library("gplots")
library("RColorBrewer")
library("genefilter")
#Read in read counts
SRR6040092 <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6040092_scaffold_10_mappedAligned_counted', header=FALSE);
colnames(SRR6040092)=c('gene', 'SRR6040092')
SRR6040093 <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6040093_scaffold_10_mappedAligned_counted', header=FALSE);
colnames(SRR6040093)=c('gene', 'SRR6040093')
SRR6040094 <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6040094_scaffold_10_mappedAligned_counted', header=FALSE);
colnames(SRR6040094)=c('gene', 'SRR6040094')
SRR6040096 <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6040096_scaffold_10_mappedAligned_counted', header=FALSE);
colnames(SRR6040096)=c('gene', 'SRR6040096')
SRR6040097 <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6040097_scaffold_10_mappedAligned_counted', header=FALSE);
colnames(SRR6040097)=c('gene', 'SRR6040097')
SRR6156066 <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6156066_scaffold_10_mappedAligned_counted', header=FALSE);
colnames(SRR6156066)=c('gene', 'SRR6156066')
SRR6156067 <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6156067_scaffold_10_mappedAligned_counted', header=FALSE);
colnames(SRR6156067)=c('gene', 'SRR6156067')
SRR6156069 <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6156069_scaffold_10_mappedAligned_counted', header=FALSE);
colnames(SRR6156069)=c('gene', 'SRR6156069')
SRR6040095_P <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6040095_scaffold_10_P_mappedAligned_counted', header=FALSE);
colnames(SRR6040095_P)=c('gene', 'SRR6040095_P')
SRR6040095_U <- read.delim('C:/Users/Mathilda Stigenberg/Dropbox/Uppsala/genomeanalysis/htseq_files/SRR6040095_scaffold_10_U_mappedAligned_counted', header=FALSE);
colnames(SRR6040095_U)=c('gene', 'SRR6040095_U')
#Merging count read data into one matrix
table_gene_counts <- merge(merge(merge(merge(merge(merge(merge(merge(merge(SRR6040092, SRR6040093, by="gene", all.x=TRUE, all.y=TRUE), SRR6040094, by="gene", all.x=TRUE, all.y=TRUE), SRR6040096, by="gene", all.x=TRUE, all.y=TRUE), SRR6040097,by="gene", all.x=TRUE, all.y=TRUE), SRR6156066, by="gene", all.x=TRUE, all.y=TRUE), SRR6156067, by="gene", all.x=TRUE, all.y=TRUE), SRR6156069, by="gene", all.x=TRUE, all.y=TRUE), SRR6040095_P, by="gene", all.x=TRUE, all.y=TRUE), SRR6040095_U, by="gene", all.x=TRUE, all.y=TRUE)
#Converting to a deseq friendly format
gene_names <- table_gene_counts$gene
sample_index <- grepl("SRR\\d+", colnames(table_gene_counts))
matrix_deseq <- as.matrix(table_gene_counts[,sample_index])
rownames(matrix_deseq) <- gene_names
#Creating matrix with tissue type
SRR <- c('SRR6040092', 'SRR6040093', 'SRR6040094', 'SRR6040096', 'SRR6040097', 'SRR6156066', 'SRR6156067', 'SRR6156069', 'SRR6040095_P', 'SRR6040095_U');
tis_type <- c('leaf', 'root', 'aril_mk', 'stem', 'aril_mk', 'aril_mon', 'aril_mon', 'aril_mon', 'aril_mk', 'aril_mk');
tissues <- data.frame(tis_type)
#Converting to a deseq friendly format
rownames(tissues) <- c(SRR)
deseq2_data <- DESeqDataSetFromMatrix(countData=matrix_deseq, colData=tissues, design= ~ tis_type)
as.data.frame(colData(deseq2_data))
dds <- DESeq(deseq2_data)
# Regularized-logarithm transform
rld <- rlog(dds)
#PCA
dev.off()
print(plotPCA(rld, intgroup=c("tis_type")))
#Heat map with gene clustering
dev.off()
topVarGenes <- head(order(rowVars(assay(rld)), decreasing=TRUE), 10)
heatmap.2(assay(rld)[topVarGenes, ], scale="row", trace="none", dendogram="column", margins=c(8,8), col=colorRampPalette(rev(brewer.pal(9, "RdBu")))(255))