Polysome Profiling in Adult Mouse Testes

Materials and reagents

1. 1.5 mL microcentrifuge tube

2. 2 mL Eppendorf tubes

3. 20 mL syringe

4. C57BL/6N male mice

5. Cycloheximide (CHX) (MCE, catalog number: HY-12320)

6. cOmplete^TM EDTA-free protease inhibitor cocktail tablets (Roche, catalog number: 04693132001)

7. SUPERase•In^TM RNase inhibitor (Invitrogen, catalog number: AM2694)

8. Triton X-100 (EMD Millipore, catalog number: 648466)

9. Sucrose, RNase- & DNase-free (AMRESCO, catalog number: 0335)

10. RNAiso Plus (Takara, catalog number: 9109)

11. Tris-HCl (1 M, pH 7.0, RNase-free) (Invitrogen, catalog number: AM9851)

12. NaCl (5 M) (Invitrogen, catalog number: AM9760G)

13. MgCl₂ (1 M) (Invitrogen, catalog number: AM9530G)

14. Chloroform (Sinopharm Chemical Reagent, catalog number: 10006818)

15. Isopropanol (Sinopharm Chemical Reagent, catalog number: 80109218)

16. Ethanol (Sinopharm Chemical Reagent, catalog number: 10009218)

17. Polysome buffer (see Recipes)

18. Lysis buffer (see Recipes)

19. 60% sucrose grading solution (see Recipes)

20. 10% sucrose grading solution (see Recipes)

21. 75% ethanol (see Recipes)

Equipment

1. Ultra-clear centrifuge tube, 14 × 89 mm (Beckman Coulter, catalog number: Z90805SCA)

2. Biocomp Gradient Maker and Fractionator (Biocomp)

3. Beckman Coulter Optima L-80 XP ultracentrifuge (with rotor SW41 Ti) (Beckman Coulter)

4. Dounce homogenizer, 2 mL (Sangon, catalog number: F519062)

5. Tweezers (Sangon Biotech, catalog number: F519022)

6. Scissor (Sangon Biotech, catalog number: F519234)

7. 100 mm cell culture dish, DNase- & RNase-free (WHB-100)

Software

We perform all sequencing data analyses under the Linux system. We use the conda program to create the environment for polysome profiling data analysis and build computational pipelines using the snakemake program with the following software:

1. Sratoolkit (https://github.com/ncbi/sra-tools/wiki/01.-Downloading-SRA-Toolkit)

2. Fastqc (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)

3. STAR (https://github.com/alexdobin/STAR)

4. Samtools (https://github.com/samtools/samtools)

5. RSeQC (https://rseqc.sourceforge.net/index.html)

6. R (https://www.r-project.org/)

7. Rsubread (https://bioconductor.org/packages/release/bioc/html/Rsubread.html)

8. DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html)

9. Xtail (https://github.com/xryanglab/xtail)

Databases

The genome sequences and gene annotation can be downloaded from the Genecode database (mm10 version 23), and tRNA and rRNA sequences can be downloaded using the UCSC Table Browser.

1. Genecode (https://www.gencodegenes.org/)

2. UCSC Genome Browser (https://genome.ucsc.edu/)

Procedure

1. Prepare mouse testicular lysate

   1. C57BL/6N male mice were sacrificed by cervical dislocation and the testes were collected after abdominal anatomy (see also Note 4).

   2. Use tweezers to peel off the tunica albuginea covering the surface of the testis.

      
      

      
      Figure 1. Schematic diagram illustrating the tunica albuginea's location and shape and remaining testicular tissues

      
      

   3. Weigh the testicular tissues and transfer them to a pre-cooled Dounce homogenizer (see Note 4).

   4. Homogenize the testicular tissues with lysis buffer (1 mL per 100 mg of testes) by 10–15 strokes.

   5. Transfer homogenate to a 1.5 mL microcentrifuge tube and incubate on ice for 15 min.

   6. Centrifuge at 13,000× g for 2 min at 4 °C.

   7. Transfer the supernatant to a fresh 1.5 mL microcentrifuge tube on ice.




2. Prepare sucrose gradient

   1. Mark the ultra-clear centrifuge tube with the marker block provided by Biocomp and place it in the MagnaBase holder.

   2. Add 10% sucrose up to the marker line with a 20 mL syringe.

   3. Then, add an equal volume of 60% sucrose at the bottom of the tube with a fresh syringe.

   4. Apply the short cap provided by Biocomp and remove the excess solution.

   5. Turn on the Biocomp Gradient Station and level the rotating platform using a round bubble level.

   6. Place the MagnaBase holder in the center of the plate.

   7. Select the pre-programmed 10%–60% (w/v) gradients and press USE.

   8. The station will beep when finished. Sucrose gradients should be used immediately.




3. Isolate polysome fractions from the total lysate

   1. Take 10% of the lysates prepared in section A to a new 1.5 mL Eppendorf tube, add 500 μL of RNAiso Plus, and mix samples well. Label this sample as TOTAL and store it at -20 °C, or extract RNA following the procedure in section D.

   2. Fill the tube prepared in section B with residual lysates.

      
      

      
      Figure 2. Loading of lysates on sucrose gradients

      
      

   3. Centrifuge the gradient at 38,000× g for 2 h at 4 °C in an SW41 Ti rotor.

   4. After ultracentrifugation, gently remove tubes from the rotor and collect 12 fractions with Piston Gradient Fractionator.

      1. Start the Gradient Station Fractionator and clean the upper valve assembly and the collection tubing.

      2. Insert the centrifuge tube into the tube holder and lock it in.

      3. Place the tube holder on the fractionator and turn it back 90° to lock it.

      4. Cut out the cap of 2 mL Eppendorf tubes and place them onto the Fraction Collector.

      5. Open the software, select the user profile, and zero the detector.

      6. After blanking, press the AIR button according to the operation manual to force air through the entire system until no liquid comes out at the end of the collection tubing.

      7. Run the fractionator and collect 12 fractions. The system will plot the absorbance in real time as the sample passes through the flow cell. Export data following software instructions.

      
      

      
      Figure 3. An absorbance graph generated by Biocomp Gradient Maker and Fractionator

      
      

   5. Transfer 700 μL of each fraction to a new 1.5 mL Eppendorf tube, add 500 μL of RNAiso Plus, and mix samples well. Label these samples as FRT#1–12 and store them at -20 °C, or extract RNA following the procedure in section D.




4. Extract polysome RNAs

   1. Vortex for 10 s and keep the homogenate at room temperature for 5 min.

   2. Centrifuge at 13,000× g for 5 min at 4 °C.

   3. Transfer the supernatant to a new centrifuge tube and add 100 μL of chloroform.

   4. Vortex for 10–15 s and keep at room temperature for 5 min.

   5. Centrifuge at 13,000× g for 15 min at 4 °C.

   6. Transfer the upper layer to a new centrifuge tube.

   7. Add one volume of isopropanol and keep it at room temperature for at least 10 min.

   8. Centrifuge at 13,000× g for 15 min at 4 °C.

   9. Remove the supernatant and wash the RNA precipitate with 1 mL of pre-cooled 75% ethanol.

   10. Centrifuge at 13,000× g for 5 min at 4 °C.

   11. Remove the supernatant and keep the RNA precipitate air dried.

   12. Dissolve the RNA precipitate with an appropriate amount of preheated nuclease-free water (see Note 5).




5. High-throughput sequencing

   1. According to the absorbance graph, combine the RNA solutions extracted from polysome fractions (usually FRT#6–12 or FRT#7–12) and label it as POLY.

   2. Subject comparable TOTAL and POLY RNAs to high-throughput sequencing with the Illumina platform (performed by Novogene Co, Ltd). Briefly, the concentration and purity of extracted RNAs were checked using Nanodrop and agarose gel electrophoresis, and integrity was further determined using Agilent 2100. Then, the libraries were constructed using NEBNext^® Ultra^TM RNA Library Prep kit for Illumina^®.

Data analysis

1. Make sure that conda is installed in your Linux system




2. Create a conda environment named PPIMT and install the snakemake and required software (sratoolkit, fastqc, STAR, samtools, rseqc, R, Rsubread, and DESeq2). Create directories named scripts and snakemake in your working directory.




3. Prepare a CSV file of sample information and the following scripts for downstream analysis (Table 1)

   
   

   Table 1. Sample list

   raw	samples	rep	cond	contrast
   SRR14139008	Hetero_polysome_rep2	rep2	Hetero_polysome	None
   SRR14139009	Hetero_polysome_rep1	rep1	Hetero_polysome	None
   SRR14139016	cKO_polysome_rep2	rep2	cKO_polysome	Hetero_polysome
   SRR14139017	cKO_polysome_rep1	rep1	cKO_polysome	Hetero_polysome
   SRR14139014	cKO_total_rep2	rep2	cKO_total	Hetero_total
   SRR14139015	cKO_total_rep1	rep1	cKO_total	Hetero_total
   SRR14139018	Hetero_total_rep2	rep2	Hetero_total	None
   SRR14139019	Hetero_total_rep1	rep1	Hetero_total	None

   
   

   1. Sample information: PP_sample_info.csv

   2. Scripts:

      1. Meta.py

         import os

         import pandas as pd

         ## input

         DATA_BASE = "../data"

         GENOME_DIR = "../../anno"

         G_INDEX_DIR = os.path.join(GENOME_DIR, "star.idx")

         r_INDEX_DIR = os.path.join(GENOME_DIR, "star.idx_rRNA")

         ANNO_DIR = GENOME_DIR

         SCRIPT_DIR = "../scripts"

         SAMPLE_INFO = "../PP_sample_info.csv"

         REF_GTF = os.path.join(GENOME_DIR, "gencode.vM23.annotation.gtf")

         REF_FA = os.path.join(GENOME_DIR, "mm10.fasta")

         rRNA_FA = os.path.join(GENOME_DIR, "mm10_rRNA.fa")

         SIZES = os.path.join(GENOME_DIR, "mm10.sizes")

         Rscript = "Rscript"

         sample_info = pd.read_csv(SAMPLE_INFO, sep=",")

         samples = sample_info['samples']

         cond = sample_info[sample_info['contrast']!="None"].cond.unique()

         ## output

         RESULTS_BASE = "../results"

         ## RNAseq

         def get_control(cond, samp_info):

         contrast = samp_info[samp_info.cond == cond].contrast.unique()

         contrast = "".join(list(contrast))

         samp = samp_info[samp_info['cond']==contrast].samples

         samp = list(samp)

         return(samp)

         def get_samp(cond, samp_info):

         samp_info = samp_info[samp_info.cond == cond]

         samp = samp_info.samples

         samp = list(samp)

         return(samp)

      2. featurecount.R

         #!/usr/bin/env Rscript

         library("Rsubread")

         args <- commandArgs(T)

         ## input

         bam <- args[1]

         gtf <- args[2]

         strandspecf <- args[3]

         name <- gsub("_uniq_sort.bam", "", basename(bam))

         base <- dirname(dirname(dirname(bam)))

         Nthread <- 8

         ## output

         DESdir <- paste(base, "featurecount", sep = "/")

         txt_feacount <- paste(DESdir, paste(name, "featureCounts.txt", sep = "_"), sep = "/")

         log_feacount <- paste(DESdir, paste(name, "featureStat.log", sep = "_"), sep = "/")

         ##run

         fc_SE <-featureCounts(

         bam,

         annot.ext = gtf,

         isGTFAnnotationFile = T,

         isPairedEnd = TRUE,

         strandSpecific = as.integer(strandspecf),

         nthreads = Nthread

         )

         write.table(fc_SE$counts,txt_feacount, quote = F,col.names = F,sep = "\t")

         write.table(fc_SE$stat,log_feacount, quote = F,col.names = F,sep = "\t")

      3. Distance.R

         #!/usr/bin/env Rscript

         ## input

         RAW_DIR <- "../results/featurecount"

         ## output

         outDir <- "../results/plot_sample_cor"

         args <- commandArgs(TRUE)

         if (length(args) == 2) {

         RAW_DIR <- args[1]

         outDir <- args[2]

         } else if (length(args) > 0) {

         stop("

         ============================================================

         Usage: 'Rscript 01_distance.R [/path/to/featurecounts_files_dir] [/out_dir]'

         Or: 'Rscript 01_distance.R' (use default parameters)

         ============================================================

         ", call.=FALSE)

         }

         library(RColorBrewer)

         library(ComplexHeatmap)

         library(circlize)

         library(dplyr)

         samples <- list.files(RAW_DIR, pattern = "_featureCounts.txt")

         t_samples <- samples[grep("total", samples)]

         p_samples <- samples[grep("polysome", samples)]

         pdf_poly <- file.path(outDir, "All_polysome_cor.pdf")

         pdf_total <- file.path(outDir, "All_total_cor.pdf")

         ## run

         plot_cor <- function(RAW_DIR, samples, pdf_cor) {

         COLS <- rev(brewer.pal(n = 5, name = "RdYlBu"))

         samp_base <- sapply(samples, function(x) strsplit(x, split = "_featureCounts")[[1]][1])

         samp_cond <- sapply(as.vector(samp_base), function(x) strsplit(x, split = "_rep")[[1]][1])

         full_name <- file.path(RAW_DIR, samples)

         tsv <- lapply(full_name, function(x) read.csv(x, sep = "\t", header = F)$V2)

         GENEID <- read.csv(full_name[1], sep = "\t", header = F)$V1

         df <- do.call(cbind, tsv)

         df <- df[apply(df, 1, sum) > 40, ]

         colnames(df) <- as.vector(samp_base)

         cor_mat <- cor(df, method = "spearman")

         print(cor_mat)

         p <- Heatmap(cor_mat, name = "Correlation", col = COLS, show_row_names = T,

         show_column_names = T, cluster_columns = T, cluster_rows = T,

         clustering_distance_rows = "spearman", clustering_distance_columns = "spearman")

         pdf(pdf_cor, width = 5, height = 4)

         print(p)

         dev.off()

         }

         plot_cor(RAW_DIR, p_samples, pdf_poly)

         plot_cor(RAW_DIR, t_samples, pdf_total)

      4. DESeq.R

         #!/usr/bin/env Rscript

         library(DESeq2)

         args <- commandArgs(T)

         ## input

         samples <- args[1]

         # samples <- gsub(" ","", samples)

         samples <- strsplit(samples, split = ",")[[1]]

         samp_cond <- strsplit(args[2], split = ",")[[1]]

         ## output

         tsv_DE <- args[3]

         ## run

         run_DESeq2 <- function(samples, samp_cond, tsv_DE) {

         samples_base <- sapply(samples, basename)

         samp_base <- sapply(samples_base, function(x) strsplit(x, split = "_featureCounts")[[1]][1])

         tsv <- lapply(samples, function(x) read.csv(x, sep = "\t", header = F)$V2)

         GENEID <- read.csv(samples[1], sep = "\t", header = F)$V1

         df <- do.call(cbind, tsv)

         rownames(df) <- GENEID

         colnames(df) <- as.vector(samp_base)

         colData <- data.frame(condition = samp_cond)

         colData$condition <- factor(samp_cond, levels = unique(samp_cond))

         dds <- DESeqDataSetFromMatrix(

         countData = df,

         colData = colData,

         design = ~condition

         )

         dds <- DESeq(dds)

         res <- results(dds)

         print(head(res))

         # write all expression gene name to a file

         all_exp <- res[res$baseMean != 0, ]

         all_exp <- as.data.frame(all_exp)

         all_exp$name <- rownames(all_exp)

         write.table(all_exp, tsv_DE, sep = "\t", quote = F, row.names = F)

         }

         run_DESeq2(samples, samp_cond, tsv_DE)

      5. polysome_volcan.R

         library(ggplot2)

         library(readr)

         ## function

         volcan_plot <- function (deg_p, name) {

         library(ggplot2)

         cutoff <- 0.58

         PValue <- 0.05

         Ylim <- 20

         pdf_Vocan <- paste(name, "volcan.pdf", sep = "_")

         deg_p <- deg_p[!is.na(deg_p$padj),]

         deg_p$padj[(-log10(deg_p$padj) > Ylim)] <- 10 ^ (-Ylim)

         ##Construct the plot object

         # with legend

         deg_p$type <- ifelse((deg_p$log2FoldChange < (-cutoff)) & (deg_p$padj < 0.05) , "Down",

         ifelse((deg_p$log2FoldChange >cutoff) & (deg_p$padj < 0.05), "Up", "NS"

         ))

         label_tab <- table(deg_p$type)

         label_stat <- paste(names(label_tab), label_tab, sep = ":")

         p <- ggplot(data = deg_p,

         aes(x = log2FoldChange,y = -log10(padj),colour = type)) +

         geom_point(size = 0.2) +

         scale_color_manual(labels=label_stat, values = c("blue", "grey", "red")) +

         geom_vline( xintercept = c(-0.58, 0.58),lty = 4,col = "grey", lwd = 0.6 ) +

         geom_hline( yintercept = -log10(0.05),lty = 4,col = "grey",lwd = 0.6) +

         xlim(c(-6, 6)) +

         ylim(c(0, Ylim)) +

         theme_bw() +

         theme(

         legend.position = c(0.8, 0.8),

         panel.grid = element_blank(),

         legend.title = element_blank())+

         labs(x = "log2(KO/WT)", y = "-log10(FDR)")

         pdf(pdf_Vocan, width = 3.5, height = 4)

         print(p)

         dev.off()

         }

         ## input

         raw_DIR <- "../results/DESeq2"

         OV <- file.path(raw_DIR, "cKO_polysome_DEG.tsv")

         ## output

         OUTDIR <- "../results/figure"

         pdf_volcan <- file.path(OUTDIR, "polysome")

         ##

         deg_p <- read_tsv(OV)

         volcan_plot(deg_p, pdf_volcan)

      6. getExp_mRNA.R

         #!/usr/bin/env Rscript

         library(readr)

         # args <- commandArgs(T)

         ## input

         RAW_DIR <- "../results/featurecount/"

         samples <- list.files(RAW_DIR, pattern = "_rpkm", full.names = T)[7:8]

         ## output

         outDir <- "../results/figure "

         tsv_mRNA <- file.path(outDir, "mRNAexp.tsv")

         ## run

         samp_base <- sapply(samples, function(x) strsplit(basename(x), split = "_rpkm")[[1]][1])

         samp_cond <- sapply(as.vector(samp_base), function(x) strsplit(x, split = "_rep")[[1]][1])

         tsv <- lapply(samples, function(x) read.csv(x, sep = "\t", header = F)$V2)

         GENEID <- read.csv(samples[1], sep = "\t", header = F)$V1

         df <- do.call(cbind, tsv)

         # df <- round(t(t(df/glt$length)))

         rownames(df) <- GENEID

         colnames(df) <- as.vector(samp_base)

         df <- df[rowSums(df) > 0, ]

         df_mean <- data.frame(

         ID = rownames(df),

         rpkm_mean = apply(df, 1, mean)

         )

         df_mean$ID <- sapply(df_mean$ID, function(x) strsplit(x, split = "\\.")[[1]][1])

         write_tsv(df_mean, tsv_mRNA)

      7. Xtail.R

         ## input

         P_DEG <- ".. /DESeq2/cKO_polysome_DEG.tsv"

         T_DEG <- ".. /DESeq2/cKO_total_DEG.tsv"

         ## output

         outDir <- "../results/figure"

         args <- commandArgs(TRUE)

         if (length(args) == 3) {

         P_DEG <- args[1]

         T_DEG <- args[2]

         outDir <- args[3]

         } else if (length(args) > 0) {

         stop("

         =============================================================

         Usage: 'Rscript 02_xtail.R [/path/to/polysome/DEG_file] [/path/to/total/DEG_file] [/out_dir]'

         Or: 'Rscript 02_xtail.R' (use default parameters)

         =============================================================

         ", call.=FALSE)

         }

         library(xtail)

         DEtsv <- file.path(outDir, "DE_TE.tsv")

         ## run

         ## function

         get_baseNmame <- function(x) {

         strsplit(basename(x), split = "_featureCounts")[[1]][1]

         }

         get_mat <- function(x) {

         count <- read.csv(x, sep = "\t", header = F)$V2

         names(count) <- read.csv(x, sep = "\t", header = F)$V1

         return(count)

         }

         ## get data matrix

         countDir <- "../results/featurecount/"

         p_samps <- list.files(countDir, full.names = T, pattern = "polysome_rep[12]_featureCounts.txt")

         t_samps <- list.files(countDir, full.names = T, pattern = "total_rep[12]_featureCounts.txt")

         p_mat <- do.call(cbind, lapply(rev(p_samps), get_mat))

         colnames(p_mat) <- sapply(rev(p_samps), get_baseNmame)

         t_mat <- do.call(cbind, lapply(t_samps, get_mat))

         colnames(t_mat) <- sapply(t_samps, get_baseNmame)

         cond <- rep(c("control", "treat"), each = 2)

         test.results <- xtail(t_mat, p_mat, cond, bins = 1000)

         res_table <- test.results$resultsTable

         res_table$ID <- row.names(res_table)

         write.table(res_table, DEtsv, sep = "\t", row.names = F, quote = F)

      8. DE_TE_scatter.R

         library(readr)

         library(ggplot2)

         library(ggrepel)

         library(clusterProfiler)

         library(org.Mm.eg.db)

         outDir <- "../results/figure"

         ## input

         DEtsv <- file.path(outDir, "DE_TE.tsv")

         P_DEG <- "../results/DESeq2/cKO_polysome_DEG.tsv"

         T_DEG <- "../results/DESeq2/cKO_total_DEG.tsv"

         ## get data matrix

         countDir <- "../results/featurecount/"

         ## output

         pdf_volcan <- file.path(outDir, "DE_TE_scatter.pdf")

         tsv_volcan <- file.path(outDir, "DE_TE_scatter.tsv")

         ## function

         get_baseNmame <- function(x) {

         strsplit(basename(x), split = "_featureCounts")[[1]][1]

         }

         get_mat <- function(x) {

         count <- read.csv(x, sep = "\t", header = F)$V2

         names(count) <- read.csv(x, sep = "\t", header = F)$V1

         return(count)

         }

         ###

         ## filter low exp genes

         test.results <- read_tsv(DEtsv)

         p_samps <- list.files(countDir, full.names = T, pattern = "polysome_rep[12]_featureCounts.txt")

         t_samps <- list.files(countDir, full.names = T, pattern = "total_rep[12]_featureCounts.txt")

         samps <- c(p_samps, t_samps)

         fpkm_mat <- do.call(cbind, lapply(samps, get_mat))

         colnames(fpkm_mat) <- sapply(samps, get_baseNmame)

         ExpGenes_bool <- apply(fpkm_mat, 1, function(x) {

         all(x > 1)

         })

         ExpGenes <- rownames(fpkm_mat)[ExpGenes_bool]

         p_DEG <- read_tsv(P_DEG)

         t_DEG <- read_tsv(T_DEG)

         useGenes <- intersect(

         p_DEG$name[!is.na(p_DEG$padj)],

         t_DEG$name[!is.na(t_DEG$padj)]

         )

         useGenes <- intersect(useGenes, ExpGenes)

         deg_p <- test.results[test.results$ID %in% useGenes, ]

         cutoff <- 0.58

         deg_p$type <- ifelse((deg_p$log2FC_TE_final < (-cutoff)) & (deg_p$pvalue.adjust < 0.05), "Down",

         ifelse((deg_p$log2FC_TE_final > cutoff) & (deg_p$pvalue.adjust < 0.05), "Up", "NS")

         )

         deg_p$Gene <- sapply(

         deg_p$ID,

         function(x) strsplit(x, split = "\\.")[[1]][1]

         )

         label_tab <- table(deg_p$type)

         label_stat <- paste(names(label_tab), label_tab, sep = ":")

         HL_lst <- c(

         "Gpd2", "Stat4",

         "Hk1", "Akap3"

         )

         eg <- bitr(deg_p$Gene, fromType = "ENSEMBL", toType = "SYMBOL", OrgDb = "org.Mm.eg.db")

         deg_p <- merge(deg_p, eg, by.x = "Gene", by.y = "ENSEMBL", all.x = T)

         p <- ggplot(deg_p, aes(x = mRNA_log2FC, y = RPF_log2FC, color = type)) +

         geom_point(size = 0.1) +

         xlab(expression("log2FC(Transcriptome)")) + # x-axis label

         ylab(expression("log2FC(Polysome)")) + # y-axis label

         geom_vline(xintercept = c(-cutoff, cutoff), lty = 4, col = "grey", lwd = 0.6) +

         geom_hline(yintercept = c(-cutoff, cutoff), lty = 4, col = "grey", lwd = 0.6) +

         scale_color_manual(labels = label_stat, values = c("blue", "grey", "red")) +

         geom_text_repel(

         data = subset(deg_p, SYMBOL %in% HL_lst),

         aes(mRNA_log2FC, RPF_log2FC, label = SYMBOL),

         box.padding = 1.0, max.overlaps = 50,

         col = "black"

         ) +

         xlim(c(-3, 3)) +

         ylim(c(-3, 3)) +

         theme_bw() +

         theme(

         legend.position = c(0.8, 0.8),

         panel.grid = element_blank(),

         legend.title = element_blank()

         )

         pdf(pdf_volcan, width = 3.5, height = 3)

         print(p)

         dev.off()

         write.table(deg_p, tsv_volcan, sep = "\t", row.names = F, quote = F)

      9. DE_TE_scatter_0.01.R

         library(readr)

         library(clusterProfiler)

         library(ggplot2)

         library(ggrepel)

         ## input

         outDir <- "../results/figure/"

         DEtsv <- file.path(outDir, "DE_TE_scatter.tsv")

         ## output

         pdf_volcan <- file.path(outDir, "DE_TE_scatter0.01.pdf")

         tsv_volcan <- file.path(outDir, "DE_TE_scatter0.01.tsv")

         ## run

         ## function

         deg_p <- read_tsv(DEtsv)

         deg_p <- deg_p[!duplicated(deg_p$Gene), ]

         cutoff <- 0.58

         pcut <- 0.01

         deg_p$type1 <- ifelse((deg_p$log2FC_TE_final < (-cutoff)) & (deg_p$pvalue.adjust < pcut), "Down",

         ifelse((deg_p$log2FC_TE_final > cutoff) & (deg_p$pvalue.adjust < pcut), "Up", "NS"))

         label_tab <- table(deg_p$type1)

         label_stat <- paste(names(label_tab), label_tab, sep = ":")

         HL_lst <- c( "Gpd2", "Stat4", "Hk1", "Akap3")

         p <- ggplot(deg_p, aes(x = mRNA_log2FC, y = RPF_log2FC, color = type1)) +

         geom_point(size = 0.1) +

         xlab(expression("log2FC(Transcriptome)")) + # x-axis label

         ylab(expression("log2FC(Polysome)")) + # y-axis label

         geom_vline(xintercept = c(-cutoff, cutoff), lty = 4, col = "grey", lwd = 0.6) +

         geom_hline(yintercept = c(-cutoff, cutoff), lty = 4, col = "grey", lwd = 0.6) +

         scale_color_manual(labels = label_stat, values = c("blue", "grey", "red")) +

         geom_text_repel(

         data = subset(deg_p, SYMBOL %in% HL_lst),

         aes(mRNA_log2FC, RPF_log2FC, label = SYMBOL),

         box.padding = 1.0, max.overlaps = 50,

         col = "black"

         ) +

         xlim(c(-3, 3)) +

         ylim(c(-3, 3)) +

         theme_bw() +

         theme(

         legend.position = c(0.8, 0.8),

         panel.grid = element_blank(),

         legend.title = element_blank()

         )

         pdf(pdf_volcan, width = 3.5, height = 3)

         print(p)

         dev.off()

         write_tsv(deg_p, tsv_volcan)




4. Prepare the references for alignment

   1. Download the directory “anno” with the required Fasta files from our polysome_profiling project, including the mouse genome (mm10) (https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/GRCm38.p6.genome.fa.gz), tRNA and rRNA sequences (generated from https://genome.ucsc.edu/cgi-bin/hgTables), and genecode v23 gtf annotation file (https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/gencode.vM23.annotation.gtf.gz).

   2. Build rRNA and tRNA index.

      STAR --runThreadN 2 --runMode genomeGenerate \

      --genomeDir r_INDEX_DIR \

      --genomeFastaFiles mm10_rRNA.fa \

      --outFileNamePrefix r_INDEX_DIR

   3. Build genome index.

      STAR --runThreadN 20 --runMode genomeGenerate \

      --genomeDir g_INDEX_DIR \

      --genomeFastaFiles mm10.fasta \

      –sjdbGTFfile anno/gencode.vM23.annotation.gtf \

      --outFileNamePrefix g_INDEX_DIR




5. Store the Fastq files of sequencing reads in a directory with the name data

   Note: If you would like to download the data presented in our original article, use the sratoolkit program and the accession number list that you can find on our PPIMT project page.

   mkdir data

   while read line

   do

   name=$line

   fastq-dump -O data --gzip $name

   done < “SRR_Acc_List.txt”




6. Check the quality of raw sequencing data using fastqc as below

   rule fastqc:

   input:

   fq1 = os.path.join(DATA_BASE, "{sample}_R1.fq.gz"),

   fq2 = os.path.join(DATA_BASE, "{sample}_R2.fq.gz"),

   output:

   fq1 = os.path.join(RESULTS_BASE, "fastqc/{sample}/{sample}_R1_fastqc.html"),

   fq2 = os.path.join(RESULTS_BASE, "fastqc/{sample}/{sample}_R2_fastqc.html"),

   shell:

   """

   fastqc -f fastq -t 2 --noextract -o {RESULTS_BASE}/fastqc/{wildcards.sample} {input.fq1} {input.fq2}

   """




7. Align the raw data to tRNA and rRNA, using the STAR program as below

   rule map_rRNA:

   input:

   R1 = os.path.join(DATA_BASE, "{sample}_R1.fq.gz"),

   R2 = os.path.join(DATA_BASE, "{sample}_R2.fq.gz"),

   rRNA_index = "%s/SAindex" % r_INDEX_DIR,

   output:

   "%s/map_rRNA/{sample}/{sample}Unmapped.out.mate1" % RESULTS_BASE,

   "%s/map_rRNA/{sample}/{sample}Unmapped.out.mate2" % RESULTS_BASE,

   threads: 14

   log:

   "%s/map_rRNA/{sample}/{sample}.log" % RESULTS_BASE,

   shell:

   """

   STAR --genomeDir {r_INDEX_DIR} --readFilesIn {input.R1} {input.R2} \

   --readFilesCommand gunzip -c --runThreadN {threads} \

   --outFileNamePrefix {RESULTS_BASE}/map_rRNA/{wildcards.sample}/{wildcards.sample} \

   --outReadsUnmapped Fastx --outFilterMatchNmin 40 --outFilterScoreMinOverLread 0 \

   --outFilterMatchNminOverLread 0 1>{log} 2>&1

   """




8. Align the unmapped reads to genome, using the STAR program as below

   rule map_genome:

   input:

   f_fastq1 = "%s/map_rRNA/{sample}/{sample}Unmapped.out.mate1" % RESULTS_BASE,

   f_fastq2 = "%s/map_rRNA/{sample}/{sample}Unmapped.out.mate2" % RESULTS_BASE,

   genome_index = "%s/SAindex" % G_INDEX_DIR,

   output:

   "%s/map_genome/{sample}/{sample}Aligned.out.bam" % RESULTS_BASE,

   threads: 14

   log:

   "%s/map_genome/{sample}/{sample}.log" % RESULTS_BASE,

   shell:

   """

   STAR --genomeDir {G_INDEX_DIR} --readFilesIn {input.f_fastq1} {input.f_fastq2} \

   --runThreadN {threads} \

   --outFileNamePrefix {RESULTS_BASE}/map_genome/{wildcards.sample}/{wildcards.sample} \

   --outSAMtype BAM Unsorted \

   --outFilterMatchNmin 40 --outFilterScoreMinOverLread 0 --outFilterMatchNminOverLread 0 \

   1>{log} 2>&1

   """




9. Filter out uniquely mapped reads to bam file and make an index using samtools program

   rule sort_index_uniqbam:

   input:

   rules.map_genome.output

   output:

   unique_bam = temp("%s/sort_index_uniqbam/{sample}/{sample}_uniq.bam" % RESULTS_BASE),

   sort_bam = "%s/sort_index_uniqbam/{sample}/{sample}_uniq_sort.bam" % RESULTS_BASE,

   bai = "%s/sort_index_uniqbam/{sample}/{sample}_uniq_sort.bam.bai" % RESULTS_BASE,

   threads: 8

   log:

   log1 = "%s/sort_index_uniqbam/{sample}/{sample}sort.log" % RESULTS_BASE,

   log2 = "%s/sort_index_uniqbam/{sample}/{sample}index.log" % RESULTS_BASE,

   shell:

   """

   {samtools} view -@ {threads} -h -bq 255 {input} > {output.unique_bam}

   {samtools} sort -@ {threads} {output.unique_bam} -o {output.sort_bam} 1>{log.log1} 2>&1

   {samtools} index {output.sort_bam} 1>{log.log2} 2>&1

   """




10. Generate bigwig files from indexed bam files, using the bam2wig.py script from rseqc package

    rule bam2bw:

    input:

    "%s/sort_index_uniqbam/{sample}/{sample}_uniq_sort.bam" % RESULTS_BASE,

    output:

    "%s/bam2bw/{sample}.bw" % RESULTS_BASE,

    threads: 1

    log:

    "%s/bam2bw/{sample}.log" % RESULTS_BASE,

    shell:

    """

    bam2wig.py -i {input} -o {RESULTS_BASE}/bam2bw/{wildcards.sample} -t 1000000000 -s {SIZES} -u \

    1>{log} 2>&1

    """




11. Count mapped reads for all annotated genes using the featureCounts program from Rsubreads package

    rule featurecount:

    input:

    bam = "%s/sort_index_uniqbam/{sample}/{sample}_uniq_sort.bam" % RESULTS_BASE,

    gtf = REF_GTF,

    output:

    "%s/featurecount/{sample}_featureCounts.txt" % RESULTS_BASE,

    "%s/featurecount/{sample}_featureStat.log" % RESULTS_BASE,

    threads: 8

    params:

    Rscript = Rscript,

    log:

    "%s/featurecount/{sample}.log" % RESULTS_BASE,

    shell:

    """

    {params.Rscript} {SCRIPT_DIR}/featurecount.R {input.bam} {input.gtf} 0 1>{log} 2>&1

    """




12. Analyze differential genes based on the reads counts of all genes in different samples with DESeq2 as below

    rule DESeq2:

    input:

    control = lambda wildcards: ["%s/featurecount/%s_featureCounts.txt" % (RESULTS_BASE, samp)

    for samp in get_control(wildcards.sample, sample_info)],

    treat = lambda wildcards: ["%s/featurecount/%s_featureCounts.txt" % (RESULTS_BASE, samp)

    for samp in get_samp(wildcards.sample, sample_info)]

    output:

    tsv = "%s/DESeq2/{sample}_DEG.tsv" % RESULTS_BASE,

    params:

    Rscript = Rscript,

    run:

    samples = ",".join(list(map(str, input.control))) + "," + ",".join(list(map(str,input.treat)))

    group = ",".join(["WT"]*2+["TREAT"]*2)

    cmd2 = " ".join([params.Rscript, SCRIPT_DIR + "/DESeq2.R", samples, group, output.tsv])

    print(cmd2)

    shell(cmd2)




13. Calculate sample correlation based on the reads counts of all genes in different samples with distance.R as below

    rule sample_cor:

    input:

    "%s/featurecount/{sample}_featureCounts.txt" % RESULTS_BASE,

    output:

    poly= "%s/plot_sample_cor/All_polysome_cor.pdf" % RESULTS_BASE,

    total= "%s/plot_sample_cor/All_total_cor.pdf" % RESULTS_BASE,

    params:

    Rscript = Rscript,

    shell:

    """

    {params.Rscript} {SCRIPT_DIR}/distance.R

    """

    
    

    
    Figure 4. Correlation analysis. The Spearman correlation coefficient among four polysome profiling experiments for Fxr1^cko and control testes with two replicates. The correlation coefficient was calculated based on reads counts of all genes between two different samples.

    
    

14. Visualize polysome profiling results with volcano plot based on the results of differential expression analysis between WT and FXR1 KO with polysome_volcan.R as below

    rule plot_polysome_volcan:

    input:

    tsv = "%s/DESeq2/{sample}_DEG.tsv" % RESULTS_BASE,

    output:

    "%s/figure/polysome_volcan.pdf" % RESULTS_BASE,

    params:

    Rscript = Rscript,

    shell:

    """

    {params.Rscript} {SCRIPT_DIR}/polysome_volcan.R

    """

    
    

    
    Figure 5. Volcano plot showing the significantly changed genes. The cutoff is FDR < 0.05 and fold-change > 1.5 between Fxr1^cko and control testes.




15. Analyze differential translational efficiency based on the reads counts of expressed genes between WT and FXR1 KO using the R package Xtail between WT and FXR1 KO as below

    rule plot_TE_scatter:

    input:

    total = "%s/DESeq2/cKO_total_DEG.tsv"" % RESULTS_BASE,

    polysome = "%s/DESeq2/cKO_poly_DEG.tsv"" % RESULTS_BASE,

    output:

    "%s/figure/DE_TE_scatter0.01.pdf" % RESULTS_BASE,

    params:

    Rscript = Rscript,

    shell:

    """

    {params.Rscript} {SCRIPT_DIR}/getExp_mRNA.R

    {params.Rscript} {SCRIPT_DIR}/Xtail.R

    {params.Rscript} {SCRIPT_DIR}/DE_TE_scatter.R

    {params.Rscript} {SCRIPT_DIR}/DE_TE_scatter_0.01.R

    """

    
    

    
    Figure 6. Differential translational efficiency analysis. Scatter plot showing mRNA level changes (x-axis) against polysome profiling changes (y-axis) between adult Fxr1^cko and control testes, with translation-deficient genes in Fxr1^cko testes shown in blue.

Notes

1. Once testes are harvested, all manipulations should be performed on ice or at 4 °C.

2. Fill the ultra-clear centrifuge tube with lysate to approximately 2–3 mm from the border to prevent the tube from collapsing during centrifugation.

3. When collecting fractions with Piston Gradient Fractionator, the collection tube needs to be flushed with ddH₂O. We recommend preheating the water to avoid clogging the tube with the sucrose solution.

4. Each sample requires two mice aged 25 days or one mouse older than 35 days. One testis of a 25-day-old mouse weighs approximately 30 mg and one testis of a 35-day-old mouse weighs approximately 40 mg.

5. The amount of water used to dissolve the RNA depends on the abundance of the RNA, or in other words, on the number of testes initially lysed. As a rule of thumb, 50–100 μL of water for INPUT yields approximately 1 μg/μL RNA, and 20–50 μL of water for FRT yields 80–800 ng/μL RNA for polysome fractions and 1–2 μg/μL RNA for RNP fractions. The volume of water used to dissolve 12 fractions should be consistent.

Recipes

1. Polysome buffer (50 mL)

   50 mM Tris-HCl (pH 7.0)

   100 mM NaCl

   5 mM MgCl₂

   Nuclease free water




2. Lysis buffer (10 mL)

   1% Triton X-100

   100 μg/mL CHX

   cOmplete^TM EDTA-free protease inhibitor cocktail

   SUPERase•In^TM RNase inhibitor

   Polysome buffer (up to 10 mL)




3. 60% (w/v) sucrose grading solution (10 mL)

   60% (w/v) sucrose

   100 μg/mL CHX

   Cocktail protease inhibitor

   RNase inhibitor

   Polysome buffer (up to 10 mL)




4. 10% (w/v) sucrose grading solution (10 mL)

   10% (w/v) sucrose

   100 μg/mL CHX

   Cocktail protease inhibitor

   RNase inhibitor

   Polysome buffer (up to 10 mL)




5. 75% (v/v) ethanol (50 mL)

   75% (v/v) ethanol

   Nuclease-free water

Acknowledgments

This protocol was originally applied in our study “LLPS of FXR1 drives spermiogenesis by activating translation of stored mRNAs” published by Science (Kang et al., 2022). We are grateful for funding by the National Key R&D Program of China (2022YFA1303300, 2022YFC2702600 and 2021YFC2700200), National Natural Science Foundation of China (3191101352, 91940305, 31830109, and 31821004), Science and Technology Commission of Shanghai Municipality (2017SHZDZX01, 19JC1410200 and 17JC1420100), Innovative research team of high-level local universities in Shanghai (SHSMU-ZDCX20210902), National Postdoctoral Program for Innovative Talent grant (BX20180331), China Postdoctoral Science Foundation grant (2018M642018), Open Fund of State Key Laboratory of Reproductive Medicine (SKLRM‐K202101), and the Foundation of Key Laboratory of Gene Engineering of the Ministry of Education. Furthermore, we would like to thank Ming-Shun Sun from the Cell Biology Core Facility in Shanghai Institute of Biochemistry and Cell Biology (SIBCB), Chinese Academy of Sciences (CAS) for technical assistance.

Competing interests

The authors declare no competing interests.

Ethics

All animal studies were approved by the Institutional Animal Care and Research Advisory Committee in SIBCB, CAS.

References

1. Chassé, H., Boulben, S., Costache, V., Cormier, P., and Morales, J. (2017). Analysis of translation using polysome profiling. Nucleic Acids Res 45(3): e15.
2. Ingolia, N. T., Brar, G. A., Rouskin, S., McGeachy, A. M. and Weissman, J. S. (2012). The ribosome profiling strategy for monitoring translation in vivo by deep sequencing of ribosome-protected mRNA fragments.Nat Protoc 7(8): 1534-1550.
3. Kang, J. Y., Wen, Z., Pan, D., Zhang, Y., Li, Q., Zhong, A., Yu, X., Wu, Y. C., Chen, Y., Zhang, X., et al. (2022). LLPS of FXR1 drives spermiogenesis by activating translation of stored mRNAs. Science 377(6607): eabj6647.
4. Sassone-Corsi, P. (2002). Unique chromatin remodeling and transcriptional regulation in spermatogenesis. Science 296(5576): 2176-2178.
5. Schwanhäusser, B., Busse, D., Li, N., Dittmar, G., Schuchhardt, J., Wolf, J., Chen, W. and Selbach, M. (2011). Global quantification of mammalian gene expression control. Nature 473(7347): 337-342.
6. Steger, K. (1999). Transcriptional and translational regulation of gene expression in haploid spermatids.Anat Embryol (Berl) 199(6): 471-487.