About adapter trimming

[nazerafauzia]

Hi Team,

I tried to trim my raw (.fq) data but failed.
My sequencing platform is DNBSEQ which uses a dual barcode adapter.
In the user manual, it mentioned that the adapter sequences should be trimmed in both 5’ ends and 3’ ends, and the sequences are:
Forward filter: AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA
Reverse filter: AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG
It didn't mention about read 1 and read 2 adapter
I input those sequences for 3' and 5' in the CutAdapt paired-end adapter trimming setup and make the same for read 1 and read 2, but failed.
Do you know what's wrong with my input?

Thanks.
Nazera

[GuyTeichman]

Hi Nazera,
Do you know if your sequencing run was single-end or paired-end? (meaning - do you get a single read per fragment, or a pair of reads per fragment)?
If your sequencing is single-end, you need to use "CutAdapt single-end adapter trimming", not paired-end trimming.

Regarding the dual barcode adapter setup:
CutAdapt has recently implemented a feature that allows you to trim a Linked Adapter, meaning a dual adapter composed of one 5' sequence and one 3' sequence. You can read more about it in the CutAdapt manual.
However, this feature is not yet supported by RNAlysis.

If you want to use RNAlysis to trim your adapters, here is what I would do:

1. Trim your reads with the 3' adapter only, discarding any untrimmed reads (by setting discard_untrimmed_reads to True).
2. Trim your already trimmed reads a second time, this time with the 5' adapter only, discarding any untrimmed reads

This way, once you are done with the two rounds of trimming, you have trimmed each direction of the reads exactly once, and only kept reads that contained both the 3' and 5' components of your linked adapter.

Please let me know if that works for you or if you have any other questions.

Good luck,
Guy.

[nazerafauzia]

Hi Guy,

Thank you for your quick reply.
Since I get a pair of reads per fragment, my sequencing run was a paired-end.
I have followed your suggestion, and it works. Thank you so much!

Now I am moving to the differential expression analysis with DESeq2 and facing an error.

Here is the log textbook:
kallisto_output_scaled_per_gene
Count matrix named 'kallisto_output_scaled_per_gene'
Warning in install.packages("BiocManager") :
'lib = "C:/Program Files/R/R-4.2.2/library"' is not writable
Error in install.packages("BiocManager") : unable to install packages
Execution halted
Installing package into 'C:/Users/nazer/OneDrive/Documents/R/win-library/4.1'
(as 'lib' is unspecified)
Error in contrib.url(repos, "source") :
trying to use CRAN without setting a mirror
Calls: install.packages -> contrib.url
Execution halted
Installing package into 'C:/Users/nazer/OneDrive/Documents/R/win-library/4.1'
(as 'lib' is unspecified)
Error in contrib.url(repos, "source") :
trying to use CRAN without setting a mirror
Calls: install.packages -> contrib.url
Execution halted
Warning in install.packages("BiocManager") :
'lib = "C:/Program Files/R/R-4.2.2/library"' is not writable
Error in install.packages("BiocManager") : unable to install packages
Execution halted

What should I do?

Thanks.
Nazera

[GuyTeichman]

Hi Nazera,
it looks like RNAlysis tried to install DESeq2, and ran into a permission error (meaning - the path where R tries to install libraries by default, "C:/Program Files/R/R-4.2.2/library"', is not writable).
If you are using Windows, you can try to either:

1. Change the properties of the directory "C:/Program Files/R/R-4.2.2/library"', so that it is no longer read-only (right click->properties->untick "read only")
2. Open RStudio as an administrator, and manually install DESeq:

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")```

Let me know if that helps. 

Good luck,
Guy. 

[GuyTeichman]

Hi Nazera,
It is possible that the differential expression analysis did not find any statistically significant genes with your given threshold (by default FDR<=0.1).
The best way to make sure would be to look at the 'pval' and 'padj' columns of your differential expression table and see what do the numbers look like.

The next things I would look at are you design matrix and comparisons - perhaps you extracted the wrong comparison, missed some samples, or mislabelled them.

Best of luck,
Guy.

[nazerafauzia]

Hi Guy,

Here is the differential expression table looks like

And the design matrix

I hope you can help me.
Thank you very much!

Anisa

[GuyTeichman]

Hi Anisa,
I don't see any obvious issue with the design matrix or results table.
However, it is noteworthy that you only have two samples per conditions (and 4 samples in total). That is a fairly small sample size, and it's therefore likely that you statistical power is fairly low and therefore DESeq2 cannot detect many differentially expressed genes.

Moreover, at least from the few genes that I can see in the results table, their expression level seems to be very low (the baseMean column). I would look at the original count matrix (using a pairplot for example, or 'describe' for descriptive statistics), to figure out whether the majority of genes are also lowly-expressed.
If so, your issue may be somewhere upstream - perhaps the alignment/feature counting went wrong and many reads were not assigned to a genomic feature, perhaps many reads were lost during adapter trimming/quality trimming, and perhaps the FASTQ files you started with had a low number of reads.

I would look into all of these direction, and try to figure out if and where things went wrong, and then re-run the rest of the analysis from that step onward.

Best of luck,
Guy.

[nazerafauzia]

Hi Guy,

Thank you so much for your reply and feedback on my data.
I got the wrong information from the company to which I sent my samples.
Actually, the data I got has been filtered out of the adapters and low-quality data, so I should just quantify it directly instead of trimming it.
Your explanation about the FASTQ files I started having a low number of reads was valid because I trimmed the data which already trimmed.

I have another thing to ask you anyway.
I need the information about the strandedness of the data to quantify using kallisto. But, something confused me when I got the answer from the company. They said the library is a stranded-specific library, and read 1 is forward, read 2 is reverse.
So, what should I choose in the stranded option since I can't choose both forward and reverse?

Thanks.
Nazera

[GuyTeichman]

Hi Nazera,
For paired-end data, 'forward' indicates the data is stranded, where the first read in the pair pseudoaligns to the forward strand of a transcript, and the second read in the pair pseudoaligns to the reverse strand of a transcript.
Therefore, I would choose 'forward' as the strandedness option for your data.

As always, it's a good practice to look at the logs in the output files and make sure that the quantification went well - meaning a respectable fraction of your reads were properly pseudoaligned to your reference transcriptome, and were then counted towards specific transcripts.

Best,
Guy.

[nazerafauzia]

Hi Guy,

Thank you so much for your answer and explanation.
I got the DEGs, and the volcano plot looks good.
But, one thing that makes me unsure about the result is among the 27423 genes from the DESeq2, I just got 123 up-regulated and 89 down-regulated genes. Is that normal?
How to make sure that a good fraction of my reads were properly pseudoaligned to my reference transcriptome and were then counted towards specific transcripts?

Thanks.
Nazera

[GuyTeichman]

Hi Nazera,
Running kallisto should generate a log file. This log file will contain statistics, such as the percentage of reads pseudoaligned/uniquely pseudoaligned. This percentage should match your expectations about the data - for example, if 50% of the reads were not aligned, something probably went wrong.

Regarding the number of DE genes - that is a plausible result, given that you have only two samples per condition, and therefore the statistical power of the DESeq2 model is fairly small. Generally speaking, a larger number of samples will allow you to more confidently characterize genes as significantly DE.

Best,
Guy.

[nazerafauzia]

Hi Guy,

I couldn't find the percentage that you mentioned.
I copied part of the text in the log and found this.

[quant] finding pseudoalignments for the reads ... done
[quant] processed 21,585,969 reads, 1,885,881 reads pseudoaligned
[quant] estimated average fragment length: 226.43
[ em] quantifying the abundances ... done
[ em] the Expectation-Maximization algorithm ran for 844 rounds

Is that mean the percentage of reads pseudoaligned in the data was low?

Thanks,
Nazera

[GuyTeichman]

Yeah, that appears like a very low alignment percentage. I think there should be additional log files in your output directory (something.log), but nevermind that.

I would guess the problem is either with your reads themselves (which you can examine using the software FastQC - it has a graphical user interface as well), or with your reference transcriptome (perhaps it is incomplete, or doesn't match your samples).

Guy