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. 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. 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. 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.
Figure 7: Heatmap after rlog transformation, showing which samples that are similar and different to each other.
Figure 8: Heatmap showing the 35 most varied genes across the samples.
Figure 9: A PCA plot showing how similar and different each sample is. One sample corresponds to a dot and the colour correspond to tissue.