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 for illumina reads and the trancriptomic reads. Below, one can observe som results from the fastQC report for pacbio as well as some results from the multiQC for illumina and the transcriptome.
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 and L75 for the corrected genome assembly.
As can be seen in table 2, the value for N50 is the contig length that produces at least half of the assembly, so using contig lengths 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. 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 GC content for the corrected genome assembly.
As one can see from figure 6, the GC content is around 33 %, and it looks like a normal distribution. It is good that there is only one peak, since that means that there is probably no contamination. The conclusion is that the quality of the assembled 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: 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. Aril_mk corresponds to aril from the musang king cultivar and aril_mon corresponds to aril from the monthong cultivar.
In figure 7, one can see that samples from leaf, root and stem are pretty similar and have low variance since these dots are close to one another in the plot. This could be due to the fact that the gene expressions similar in these tissues and that is why these dots vary only a little. Similarly, if one look at the three dots close to each other that corresponds to two aril from monthong and one aril from the musang king, one can conclude that the expression for those are quite similar. However, these are far away from the rest of the dots, which may be due to the fact that these dots express more than for example the tissues leaf, stem and root. Furthermore, two dots are also close to each other, and these are arils from musang king. These probably also have similar expression level. There is also one dot corresponding to aril from musang king that is far away from all the other dots, which propably means that this sample has a high variance compared to the other samples. Similarly, there is also one dot corresponding to aril from monthong that is far away from the other, and the conclusion is that this sample also differ a lot in expression level in comparsion with the others. Since all the arils are pretty far away from one another, the conclusion is that maybe taste and smell differ from fruit to fruit.
Figure 8: Heatmap showing the 10 most varied genes across the samples.
In figure 8, one can see a heatmap of the genes with the highest variance. The gene g2884.t1.exon4 has a has a dark red colour in SRR6040095_U, which corresponds to the tissue aril from the musang king cultivar. It means that gene has a high variance in expression for that tissue compared to the others. Another gene that vary much in expression is g2734.t1.exon1 in SRR6156069, which corresponds to the tissue aril from the monthong cultivar. Two genes that vary sligthly less in express but really much anyway are g3157.t1.exon3 and g302.t1.exon1. G3157.t1.exon3 has a varied expression level in SRR6156066, which correspond to tissue aril from the monthong cultivar. G302.t1.exon1 has a varied expression level 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).