Analyses
Table 1: The metadata, information about the data.
In table 1, the metadata can be observed. All the SRR numbers per tissue, type of sample, type of tissue and the sequencing technique can be seen.
The raw data obtained were quality checked by using the software FastQC. The commands used for this step can be seen in the code page section fastaqc. Thereafter, a scaffold was chosen to continue the other analyses with. The scaffold 10 was chosen and a script was done using the program multiQC (the script can be seen in the code page section multiqc) so that the quality reports for scaffold 10 were gathered into a summary report instead. The generated summary reports from illumina, pacbio and the trancriptome were downloaded and below can some of the results from the multiQC be seen.
Illumina quality report
Figure 1: A graph of the mean quality scores of the illumina reads.
As one can see in figure 1, the quality of the illumina reads are quite high, since the curves are located in the green area. These might have a high quality due to the fact that illumina reads have a low error rate. It is important that these reads have a high quality since the reads are later going to be used for correcting the genome assembly by mapping the reads to the genome assembly and use pilon. Since the quality is good, the correcting of the genome assembly will probably be good as well.
Pacbio quality report
Figure 2: A graph of the mean quality score of the pacbio read.
From figure 2, one can conclude that the quality of the pacbio read across the base pairs is pretty bad since the curve is located in the red area. This can be due to the fact that pacbio read has a higher error rate than the illumina reads.
Transcriptomic quality report
Trimmed:
Figure 3: A graph of the mean quality scores of the trimmed transcriptomic reads.
The trancriptomic reads have high quality across the base pairs since the curves are located in the green area, as one can see in figure 3 above. This high quality is important for the analysis step when the reads are supposed to be mapped against the genome assembly. This will probably result in a good mapping since the quality is high!
Untrimmed:
Figure 4: A graph of the mean quality scores of the untrimmed transcriptomic reads.
As one can see in the figure 4 above, the quality of the untrimmed transcriptomic reads across the base pairs are quite high. The quality is good since the curves are located in the green area.
Trimmed data that previously was untrimmed:
Figure 5: A graph of the mean quality scores of the newly trimmed transcriptomic reads.
In this figure, figure 5, one can see the newly trimmed transcriptomic reads that have been trimmed with trimmomatic. The quality of these trimmed transcriptomic reads are quite high since the curves are located in the green area. If comparing with figure 4, the quality is quite similar for 2 of the reads. However, this graph also have 2 additional reads which have a lower quality than the other 2 in this graph. This may be due to the fact that the reads are now cropped compared to before, since they are now trimmed.
The pacbio scaffold 10 data was assembled using the program canu. A script was made and can be seen in the code page section canu and was run with the command sbatch.
The untrimmed transcriptomic data was trimmed using the program trimmomatic. A script was made and can be seen in the code page section trimmomatic and was run with the command sbatch. This was done in order to crop and trim the reads, as well as remove adapters. This lead to better results for the future analyses, since otherwise adapters can cause some problems.
Mapping of the illumina reads (scaffold 10 1P and scaffold 10 2P) against the assembled pacbio scaffold 10 data was done with the program BWA so that the assembly could be corrected with the illumina reads. The code can be seen in the code page section mapping with bwa and the script was run with the command sbatch.
Improvement of the assembly was made by using the program pilon that required the assembly as imput as well as the aligned bam file from the previous step. A script was made and can be seen in the code page section pilon.
Quality assessment was made with QUAST and the output from the previous analysis (pilon) was used as input. The script can be seen in the code page section quast.
Table 2: A table showing the values for N50, N75, L50, L75 and the GC content for the corrected genome assembly.
As can be seen in table 2, the value for N50 is the contig lenght that produces at least half of the assembly, so using contig length of at least 276 902 bp produces 50 % of the assembly. N75 is similar to N50 but accounts for 75 % of the assembly instead, so using contig lengths of at least 29 206 bp produces 75 % of the assembly. The value for L50 is 20, which means that 20 contigs are needed to produce at least half of the assembly. The value for L75 is 114, so 114 contigs are needed to at least produce 75 % of the assembly. The GC content is 31,44 %, so there is 31,44 % of G and C bases in the assembled genome. None of these values seems totally out of the blue, so the assumption is that the quality of the assembled genome is quite good.
Figure 6: The figure is visualizing the contig sizes for the corrected genome assembly.
As one can see from figure 6, the contig size increases. There are both short and long contigs in the corrected assembled genome and it looks like there is a normal distribution of short and long reads. The plot is cumulative and there is no sudden peak, so one can conclude that the quality of the corrected genome is good.
The program repeatMasker was used on the fasta file obtained from the pilon analysis. This program was used in order to improve the prediction accuracy by masking the repeats in the genome. The type of masking used was hardmasking, which means that the letter N replaced repetitive regions. A script was made and can be seen in the code page section repeatmasker.
The program star was used in order to map the trimmed RNA-seq reads to the assembled masked genome. A script was made and can be seen in the code page section star.
BRAKER was used to do a structural annotation of the genome. The masked genome obtained from the repeatMasker analysis and all bam files from the previous step when star was used were inputs when doing this analysis. A script was made and can be seen in the code page section braker.
The program eggNOGmapper was used to do a functional annotation of the assembled masked genome obtained from the repeatMasker analysis. A script was made and can be seen in the code page section eggnogmapper.
HTseq was used to count how many genes each transcriptomic read mapped to. As a result from this analysis, one can determine the expression level of specific genes. The BAM-files from the program star was used as input as well as a gff3 file from BRAKER. A script was made and can be seen in the code page section htseq.
A script was done in R where the packages DESeq2, gplot, genefilter and RcolorBrewer were used in order to do some gene expression analyses with the output from HTseq. The output from HTseq was the amount of gene counts each transcriptomic read had. The expression analyses were visualized with one heatmap and one PCA plot. The script can be seen in the code page section deseq2.
Figure 7: Heatmap showing the 10 most varied genes across the samples.
In figure 7, one can see a heatmap of the genes with the highest variance. The gene g2884.t1.exon4 express really much since it has a dark red colour in SRR6040095_U, which corresponds to the tissue aril from the musang king cultivar. Another gene that express much is g2734.t1.exon1 in SRR6156069, which corresponds to the tissue aril from the monthong cultivar. Two genes that express slightly less but really much anyway are g3157.t1.exon3 and g302.t1.exon1. G3157.t1.exon3 express much in SRR6156066, which correspond to tissue aril from the monthong cultivar. G302.t1.exon 1 express much in SRR6156067, which corresponds to tissue aril from the monthong cultivar. When looking in the augustus.hints.gff3 file from BRAKER, one can see that the gene g2884.t1.exon4 corresponds to tig00000563_pilon, gene g2734.t1.exon1 corresponds to tig00000073_pilon, gene g3157.t1.exon3 corresponds to tig00000603_pilon and gene g302.t1.exon1 corresponds to tig00000005_pilon. In order to find out the function these genes have, one can look in the functional_annotation.emapper.annotations file from eggNOGmapper and search on the pilon names with the grep command (etc. grep "tig00000005_pilon" functional_annotation.emapper.annotations).
Figure 8: A PCA plot showing how similar and different each sample are. One sample corresponds to a dot and the colours correspond to a tissue type.