Example: differential expression analysis

Example: RNA-Seq differential expression patterns determination

To determine statistically significant k-mers of the African Turquoise Killifish liver tissue for sex, the user may prepare an input samples.txt file of samples:

SRR22013784 SRR22013784_1_val_1.fq.gz SRR22013784_2_val_2.fq.gz
SRR22013785 SRR22013785_1_val_1.fq.gz SRR22013785_2_val_2.fq.gz
SRR22013786 SRR22013786_1_val_1.fq.gz SRR22013786_2_val_2.fq.gz
SRR22013780 SRR22013780_1_val_1.fq.gz SRR22013780_2_val_2.fq.gz
SRR22013769 SRR22013769_1_val_1.fq.gz SRR22013769_2_val_2.fq.gz
SRR22013763 SRR22013763_1_val_1.fq.gz SRR22013763_2_val_2.fq.gz
SRR22013779 SRR22013779_1_val_1.fq.gz SRR22013779_2_val_2.fq.gz
SRR22013793 SRR22013793_1_val_1.fq.gz SRR22013793_2_val_2.fq.gz
SRR22013782 SRR22013782_1_val_1.fq.gz SRR22013782_2_val_2.fq.gz
SRR22013774 SRR22013774_1_val_1.fq.gz SRR22013774_2_val_2.fq.gz
SRR22013762 SRR22013762_1_val_1.fq.gz SRR22013762_2_val_2.fq.gz
SRR22013770 SRR22013770_1_val_1.fq.gz SRR22013770_2_val_2.fq.gz
SRR22013775 SRR22013775_1_val_1.fq.gz SRR22013775_2_val_2.fq.gz
SRR22013765 SRR22013765_1_val_1.fq.gz SRR22013765_2_val_2.fq.gz
SRR22013787 SRR22013787_1_val_1.fq.gz SRR22013787_2_val_2.fq.gz
SRR22013789 SRR22013789_1_val_1.fq.gz SRR22013789_2_val_2.fq.gz

prepare design/phenotype file liver_sex.txt:

male
male
female
female
female
female
female
female
male
male
male
female
male
male
female
male

and then perform dimensionality reduction with the following command line:

./mkmc                                           (1)
 -k 30 -m 16 \                                   (2)
--thr-rat 0.5 --thr 1 \                          (3)
--flt GCF_001465895.1_Nfu_20140520_genomic.fna \ (4)
 -t 32 \                                         (5)
--cs 1000000000 \                                (6)
-n freq \                                        (7)
--diff wrs -c liver_sex.txt \                    (8)
--pval-corr bh                                   (9)
--max-corrected-pval 0.05 \                      (10)
-- samples.txt GSE216369_liver_trimmed tmp       (11)

The consecutive lines mean as follows:

1. Run MKMC.
2. Count 30-mers; k-mer counting should not use more than 16 GB of memory; by default, samples data is stored in FASTQ files (see Building a matrix for non-FASTQ input files).
3. Exclude k-mers from the matrix which ones has no counts in at least half of the samples (see Filtering k-mers along with building the matrix).
4. Keep only k-mers which are present in GCF_001465895.1_Nfu_20140520_genomic.fna file (see Filtering k-mers along with building the matrix). Here we pass the file with the whole genome to consider k-mers only present there.
5. Use 32 threads.
6. Set the maximal number of count to $10^9$ instead of $65535$.
7. Normalize counts with frequency count method, but do not save them to the disk (rather use them as Wilcoxon-rank sum (WRS) input, see Normalizing counts).
8. Perform differential k-mers expression analysis with WRS. The design/phenotype is stored in liver_sex.txt file (see Performing differential k-mers analysis).
9. Correct WRS p-values with Benjamini-Hochberg method (see Correcting p-values).
10. Treat k-mers with p-value no higher than 0.05 as statistically significant (see Correcting p-values).
11. Pass to MKMC the file with samples, an output file name template (GSE216369_liver_trimmed) and a temporary directory (tmp). Note the space after the two dashes.

The output GSE216369_liver_trimmed_wrs_cor_significant.fa file contains statistically significant k-mers. It may be utilized to extract significant genes by mapping k-mers (e.g. with BWA), summarizing the mappings (featureCounts), and extracting the genes to a file GSE216369_liver_trimmed.significant:

bwa mem -t 32 -O 1000 -L 1000 -B 10000 -k 30 GCF_001465895.1_Nfu_20140520_genomic.fna GSE216369_liver_trimmed_wrs_cor_significant.fa > GSE216369_liver_trimmed_aligned.sam
featureCounts -T 32 -a ref_Nfu_20140520_top_level_filtered.gtf -g gene_name -o GSE216369_liver_trimmed.featureCounts GSE216369_liver_trimmed_aligned.sam
awk 'NR>2 { if ($7 > 0) print $1; }' GSE216369_liver_trimmed.featureCounts | sort > GSE216369_liver_trimmed.significant

To trim the input reads (as in the example), the user may use e.g. trim_galore:

./trim_galore -j 8 --fastqc --fastqc_args "--outdir fastqc/SRR22013784" --gzip --output_dir . --paired SRR22013784_1.fastq.gz SRR22013784_2.fastq.gz