引言

❝
相关介绍 Science 丨“千里姻缘一线牵”——蛋白质组装的秘密。德国海德堡大学的 Günter Kramer 和 Bernd Bukau 教授在 2021 年 Science 上提出了 共翻译组装（co-co assembly） 的核糖体组装机制，该研究通过蛋白质组分析，发现共翻译-组装（co-co assembly）是新生蛋白亚基组装成蛋白二聚体的通用且高效的策略。

通过 DISP（Disome selective profiling） 在全蛋白组揭示新生肽的互相作用：

分析流程：

❝
作者发现了一种先前未被认识的蛋白质复合物组装机制，其中相邻核糖体翻译的 nascent 蛋白质之间的相互作用（co-co 组装）驱动了同源复合物的形成。使用一种称为 disome 选择性分析（DiSP）的全蛋白质组筛查方法，作者在不同细胞亚细胞结构中鉴定了数千个 undergo co-co 组装的候选蛋白质，主要形成同源复合物。co-co 组装过程由五种保守的结构 motif 介导，其中卷曲螺旋是最常见的。作者能够在大肠杆菌中重建核层蛋白的 co-co 组装，证明这一过程可以独立于专门的真核细胞组装因子发生，主要依赖 nascent 蛋白质 N 端的二聚化倾向。作者提出，多聚体核糖体为大多数 co-co 组装相互作用提供了平台，提高了同源复合物形成的效率和准确性，使细胞能够独立进化功能多样的同源蛋白质复合物。
数据下载

GSE151959：

下载方式多样，自己选择即可。
去接头

❝
注意这里不同测序平台的接头不一样，所以需要分开进行处理。

trim adpator

for i in Dis_Ecoli_1B-Cherry_Rep1 Dis_Ecoli_1B_Cherry_Rep1 Dis_Ecoli_1B-Cherry_Rep2 Dis_Ecoli_1B_Cherry_Rep2 Dis_Ecoli_DCTN1_Rep1 Dis_Ecoli_LMNC_Rep1 Dis_Ecoli_LMNC_Rep2 Dis_HEK_no_PK Dis_HEK_no_puro_Rep1 Dis_HEK_no_puro_Rep2 Dis_HEK_PK_high Dis_HEK_PK_low Dis_HEK_PK_mid Dis_HEK_PK_veryhigh Dis_HEK_puro_Rep1 Dis_HEK_puro_Rep2 Dis_HEK_Rep1 Dis_HEK_Rep2 Mono_Ecoli_1B-Cherry_Rep1 Mono_Ecoli_1B_Cherry_Rep1 Mono_Ecoli_1B-Cherry_Rep2 Mono_Ecoli_1B_Cherry_Rep2 Mono_Ecoli_DCTN1_Rep1 Mono_Ecoli_LMNC_Rep1 Mono_Ecoli_LMNC_Rep2 Mono_HEK_no_PK Mono_HEK_no_puro_Rep1 Mono_HEK_no_puro_Rep2 Mono_HEK_PK_high Mono_HEK_PK_low Mono_HEK_PK_mid Mono_HEK_PK_veryhigh Mono_HEK_puro_Rep1 Mono_HEK_puro_Rep2 Mono_HEK_Rep1 Mono_HEK_Rep2
do
cutadapt -q 20 -m 23 -j 20 –discard-untrimmed -O 6 \
-a ATCGTAGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
-o ./${i}_trimed.fastq.gz \
./${i}_umi.fastq.gz
done

code for U2OS cell

for i in Dis_U2OS_Rep1 Dis_U2OS_Rep2 Mono_U2OS_Rep1 Mono_U2OS_Rep2
do
cutadapt -q 20 -m 23 -j 20 –discard-untrimmed -O 6 \
-a CTGTAGGCACCATCAATTCGTATGCCGTCTTCTGCTTG \
-o ./${i}_trimed.fastq.gz \
./${i}_umi.fastq.gz
done
处理 UMI 序列

❝
reads 两端含有 UMI 序列，进行去除然后添加到名称里面：

add UMI

for i in Dis_Ecoli_1B-Cherry_Rep1 Dis_Ecoli_1B_Cherry_Rep1 Dis_Ecoli_1B-Cherry_Rep2 Dis_Ecoli_1B_Cherry_Rep2 Dis_Ecoli_DCTN1_Rep1 Dis_Ecoli_LMNC_Rep1 Dis_Ecoli_LMNC_Rep2 Dis_HEK_no_PK Dis_HEK_no_puro_Rep1 Dis_HEK_no_puro_Rep2 Dis_HEK_PK_high Dis_HEK_PK_low Dis_HEK_PK_mid Dis_HEK_PK_veryhigh Dis_HEK_puro_Rep1 Dis_HEK_puro_Rep2 Dis_HEK_Rep1 Dis_HEK_Rep2 Dis_U2OS_Rep1 Dis_U2OS_Rep2 Mono_Ecoli_1B-Cherry_Rep1 Mono_Ecoli_1B_Cherry_Rep1 Mono_Ecoli_1B-Cherry_Rep2 Mono_Ecoli_1B_Cherry_Rep2 Mono_Ecoli_DCTN1_Rep1 Mono_Ecoli_LMNC_Rep1 Mono_Ecoli_LMNC_Rep2 Mono_HEK_no_PK Mono_HEK_no_puro_Rep1 Mono_HEK_no_puro_Rep2 Mono_HEK_PK_high Mono_HEK_PK_low Mono_HEK_PK_mid Mono_HEK_PK_veryhigh Mono_HEK_puro_Rep1 Mono_HEK_puro_Rep2 Mono_HEK_Rep1 Mono_HEK_Rep2 Mono_U2OS_Rep1 Mono_U2OS_Rep2
do
./abc7151_Script_1.jl ../1.raw-data/${i}_trimed.fastq.gz ./${i}_umi.fastq.gz –umi3 5 –umi5 2
done
remove rRNA

❝
下载 人 和 大肠杆菌 的 rRNA 然后建立索引，使用 bowtie2 进行去除：

remove Ecoli rRNA

for i in Dis_Ecoli_1B-Cherry_Rep1 Dis_Ecoli_1B_Cherry_Rep1 Dis_Ecoli_1B-Cherry_Rep2 Dis_Ecoli_1B_Cherry_Rep2 Dis_Ecoli_DCTN1_Rep1 Dis_Ecoli_LMNC_Rep1 Dis_Ecoli_LMNC_Rep2 Mono_Ecoli_1B-Cherry_Rep1 Mono_Ecoli_1B_Cherry_Rep1 Mono_Ecoli_1B-Cherry_Rep2 Mono_Ecoli_1B_Cherry_Rep2 Mono_Ecoli_DCTN1_Rep1 Mono_Ecoli_LMNC_Rep1 Mono_Ecoli_LMNC_Rep2
do
bowtie2 -p 20 -L 13 -x ../../index-data/Ecoli-rRNA-index/Ecoli-rRNA \
–un-gz ./${i}.rmrRNA.fq.gz \
-U ../2.trim-data/${i}_umi.fastq.gz \
-S ./null
done

remove human rRNA

for i in Dis_HEK_no_PK Dis_HEK_no_puro_Rep1 Dis_HEK_no_puro_Rep2 Dis_HEK_PK_high Dis_HEK_PK_low Dis_HEK_PK_mid Dis_HEK_PK_veryhigh Dis_HEK_puro_Rep1 Dis_HEK_puro_Rep2 Dis_HEK_Rep1 Dis_HEK_Rep2 Dis_U2OS_Rep1 Dis_U2OS_Rep2 Mono_HEK_no_PK Mono_HEK_no_puro_Rep1 Mono_HEK_no_puro_Rep2 Mono_HEK_PK_high Mono_HEK_PK_low Mono_HEK_PK_mid Mono_HEK_PK_veryhigh Mono_HEK_puro_Rep1 Mono_HEK_puro_Rep2 Mono_HEK_Rep1 Mono_HEK_Rep2 Mono_U2OS_Rep1 Mono_U2OS_Rep2
do
bowtie2 -p 20 -L 13 -x ../../index-data/human-rRNA-index/human-rRNA \
–un-gz ./${i}.rmrRNA.fq.gz \
-U ../2.trim-data/${i}_umi.fastq.gz \
-S ./null
done
比对到基因组

❝
分别下载 人 和 大肠杆菌 的基因组和注释文件建立索引，使用 STAR 软件进行比对。大肠杆菌的网址：http://bacteria.ensembl.org/Escherichia_coli_gca_900608175/Info/Index/

这里还需要使用 gffread 将 gff3 转为 gtf 格式来给 star 建索引。

mapping for human

index

mkdir GRCh38-star-genome-index

STAR –runThreadN 20 \
–runMode genomeGenerate \
–genomeDir GRCh38-star-genome-index/ \
–genomeFastaFiles ./Homo_sapiens.GRCh38.dna.primary_assembly.fa \
–sjdbGTFfile Homo_sapiens.GRCh38.112.gtf

mapping for Ecoli

index

gff3 to gtf

gffread Escherichia_coli_gca_900608175.ecoli025.59.gff3 -T -o Escherichia_coli_gca_900608175.ecoli025.59.gtf

mkdir Ecoli-star-genome-index

STAR –runThreadN 20 \
–runMode genomeGenerate \
–genomeDir Ecoli-star-genome-index/ \
–genomeFastaFiles ./Escherichia_coli_gca_900608175.ecoli025_.dna.toplevel.fa \
–sjdbGTFfile Escherichia_coli_gca_900608175.ecoli025.59.gtf
然后分别比对：

mapping

for i in Dis_HEK_no_PK Dis_HEK_no_puro_Rep1 Dis_HEK_no_puro_Rep2 Dis_HEK_PK_high Dis_HEK_PK_low Dis_HEK_PK_mid Dis_HEK_PK_veryhigh Dis_HEK_puro_Rep1 Dis_HEK_puro_Rep2 Dis_HEK_Rep1 Dis_HEK_Rep2 Dis_U2OS_Rep1 Dis_U2OS_Rep2 Mono_HEK_no_PK Mono_HEK_no_puro_Rep1 Mono_HEK_no_puro_Rep2 Mono_HEK_PK_high Mono_HEK_PK_low Mono_HEK_PK_mid Mono_HEK_PK_veryhigh Mono_HEK_puro_Rep1 Mono_HEK_puro_Rep2 Mono_HEK_Rep1 Mono_HEK_Rep2 Mono_U2OS_Rep1 Mono_U2OS_Rep2
do
STAR –runThreadN 20 \
–outFilterType Normal \
–alignIntronMin 5 \
–outFilterMismatchNmax 1 \
–outFilterMultimapNmax 1 \
–genomeDir ../../index-data/GRCh38-star-genome-index \
–readFilesCommand gunzip -c \
–readFilesIn ../3.rmrRNA-data/${i}.rmrRNA.fq.gz \
–outFileNamePrefix ./${i} \
–outSAMtype BAM SortedByCoordinate \
–quantMode GeneCounts \
–outSAMattributes All XS
done

mapping

for i in Dis_Ecoli_1B-Cherry_Rep1 Dis_Ecoli_1B_Cherry_Rep1 Dis_Ecoli_1B-Cherry_Rep2 Dis_Ecoli_1B_Cherry_Rep2 Dis_Ecoli_DCTN1_Rep1 Dis_Ecoli_LMNC_Rep1 Dis_Ecoli_LMNC_Rep2 Mono_Ecoli_1B-Cherry_Rep1 Mono_Ecoli_1B_Cherry_Rep1 Mono_Ecoli_1B-Cherry_Rep2 Mono_Ecoli_1B_Cherry_Rep2 Mono_Ecoli_DCTN1_Rep1 Mono_Ecoli_LMNC_Rep1 Mono_Ecoli_LMNC_Rep2
do
STAR –runThreadN 20 \
–outFilterType Normal \
–alignIntronMin 5 \
–outFilterMismatchNmax 1 \
–outFilterMultimapNmax 1 \
–genomeDir ../../index-data/Ecoli-star-genome-index \
–readFilesCommand gunzip -c \
–readFilesIn ../3.rmrRNA-data/${i}.rmrRNA.fq.gz \
–outFileNamePrefix ./${i} \
–outSAMtype BAM SortedByCoordinate \
–quantMode GeneCounts \
–outSAMattributes All XS
done
assign E/P/A site

❝
然后给比对的文件分配 EPA 位点，并输出 HDF5 文件为后续分析。因为代码本身是有时效性的，不同版本的更新迭代，会使很多代码使用不了，比如这篇文章 2021 年的，使用他提供的 julia 脚本也许会发生很多意想不到的错误。

debug 过程：
脚本前面这些导入的都需要确定安装好了：

单独运行脚本遇到报错，原因是脚本并没有导入 XAM：

解决以上报错后，单独运行脚本遇到报错，找不到 DNAKmer：

解决以上报错后，单独运行脚本遇到报错，找不到 GFF3：

所以我们需要安装这些包，然后在脚本里进行导入：

!/usr/bin/env julia

Base.CoreLogging.disable_logging(Base.CoreLogging.Debug)

using BioAlignments
using BioSequences
using GenomicFeatures
using DataStructures

using XAM
using Kmers
using GFF3

using CodecZlib
using HDF5

using Gadfly
import Cairo
import Compose

using ArgParse
解决以上报错后，单独运行脚本遇到报错，找不到 DNASequences：

根据报错信息找到相关位置：

去看看 BioSequences 说明文档怎么构造序列的，发现并不是像上面这样 DNASequence 使用的：

我们需要修改成这样：

umiseq = DNASequence(umistr[(findlast(“__“, umistr).stop + 1):end])

tmp_seq = umistr[(findlast(“__“, umistr).stop + 1):end]
umiseq = LongDNA{4}(string(tmp_seq))
继续测试，报错 g_create：

寻找报错位置：

g_create 应该是创建 group dataset，HDF5 里则应该是使用 create_group 来完成：

修改成这样：

g = Dict((i, g_create(f, string(rsite))) for (i, f) ∈ out_counts)

g = Dict((i, create_group(f, string(rsite))) for (i, f) ∈ out_counts)
继续测试，报错找不到 Name：

寻找相关位置，前面出现的地方都是条件判断，前者为真就取前面的，标红这个地方是直接取，对于取不到的情况也许就会报错：

我们打开 GFF3 文件找到 gene 的地方，有一些基因就是没有 Name 的，所以就出现了上面取不到值报错的情况，我们添加测试的打印代码看看到哪个基因出现报错：

然后找到这个基因，往下一个基因就是有问题的，可以看到也符合我们的预期：

使用 try catch 来捕捉这些错误，对于不满足条件的情况我们直接取 gene id 就行了，然后顺便打印出哪些基因不满足条件：
elseif GFF3.featuretype(record) == “gene”
att = Dict(GFF3.attributes(record))
id = haskey(att, “ID”) ? att[“ID”][1] : att[“Name”][1]
chrom = GFF3.seqid(record)

# println(att["Name"][1])
# gene = Gene(att["Name"][1], id, chrom)

gene = Gene("", "", "")

try
    gene = Gene(att["Name"][1], id, chrom)
catch e
    @info id*" has no name! Using gene id instead."
finally
    gene = Gene(att["ID"][1], id, chrom)
end

for k ∈ keys(GeneAttributes)
        if haskey(att, k)
            gene.attributes[k] = att[k][1]
        end
    end
genes[id] = gene

end
继续测试，报错提示找不到 d_create：

和之前类似的，是 HDF5 创建 dataset 的问题：

需要修改成这样：

dset = d_create(h5out, gene.name, arr)[1]

dset = create_dataset(h5out, gene.name, arr)[1]
修改完后就可以正常输出文档了：

最后进行分析：

wget https://ftp.ensembl.org/pub/release-112/gff3/homo_sapiens/Homo_sapiens.GRCh38.112.chr.gff3.gz

gunzip Homo_sapiens.GRCh38.112.chr.gff3.gz

assignment for human

for i in Dis_HEK_no_PK Dis_HEK_no_puro_Rep1 Dis_HEK_no_puro_Rep2 Dis_HEK_PK_high Dis_HEK_PK_low Dis_HEK_PK_mid Dis_HEK_PK_veryhigh Dis_HEK_puro_Rep1 Dis_HEK_puro_Rep2 Dis_HEK_Rep1 Dis_HEK_Rep2 Dis_U2OS_Rep1 Dis_U2OS_Rep2 Mono_HEK_no_PK Mono_HEK_no_puro_Rep1 Mono_HEK_no_puro_Rep2 Mono_HEK_PK_high Mono_HEK_PK_low Mono_HEK_PK_mid Mono_HEK_PK_veryhigh Mono_HEK_puro_Rep1 Mono_HEK_puro_Rep2 Mono_HEK_Rep1 Mono_HEK_Rep2 Mono_U2OS_Rep1 Mono_U2OS_Rep2
do
./abc7151_Script_2.jl -c 1 -g Homo_sapiens.GRCh38.112.chr.gff3 -u -o ${i} ../4.map-data/${i}Aligned.sortedByCoord.out.bam
done

assignment for Ecoli

for i in Dis_Ecoli_1B-Cherry_Rep1 Dis_Ecoli_1B_Cherry_Rep1 Dis_Ecoli_1B-Cherry_Rep2 Dis_Ecoli_1B_Cherry_Rep2 Dis_Ecoli_DCTN1_Rep1 Dis_Ecoli_LMNC_Rep1 Dis_Ecoli_LMNC_Rep2 Mono_Ecoli_1B-Cherry_Rep1 Mono_Ecoli_1B_Cherry_Rep1 Mono_Ecoli_1B-Cherry_Rep2 Mono_Ecoli_1B_Cherry_Rep2 Mono_Ecoli_DCTN1_Rep1 Mono_Ecoli_LMNC_Rep1 Mono_Ecoli_LMNC_Rep2
do
./abc7151_Script_2.jl -c 1 -g ../../index-data/Escherichia_coli_gca_900608175.ecoli025.59.gff3 -u -o ${i} ../4.map-data/${i}Aligned.sortedByCoord.out.bam
done
最后每个文件夹会产生这些文件，作者是拿 P-site 文件进行后续分析的：

作者也上传了 HDF5 文件，后续我们拿这个来分析：

metagene analysis

❝
作者使用了 position-wise 的方法来计算 density/enrichment。假定 counts 的分布符合泊松分布，且 disome 和 monosome 是随机独立的。然后使用 Agresti-Coull method 方法来计算置信区间。为了减小一些虚假峰的影响，使用 15 个 codon 的滑动窗口。

❝
这里作者使用 RiboSeqTools 来做后续的分析，应该是自己写的，但是缺少使用文档，我三年前问的，让作者写个使用文档，到现在也没有回复，没办法，我们只能花时间取把每个代码都测试看看：

❝
使用 load_serp 函数将 hdf5 数据导入进来，整理为一个 S3 对象，后续大部分函数都是基于这个对象操作的：
devtools::install_github(“ilia-kats/RiboSeqTools”)

library(RiboSeqTools)

make object

serp <- load_serp(HEK =list(ip = c(“GSE151959_RAW/GSM4594341_Dis_HEK_Rep1.h5”,
“GSE151959_RAW/GSM4594343_Dis_HEK_Rep2.h5”),
tt = c(“GSE151959_RAW/GSM4594340_Mono_HEK_Rep1.h5”,
“GSE151959_RAW/GSM4594342_Mono_HEK_Rep2.h5”)),
# normalize = TRUE,
bin = ‘byaa’)

metagene_profiles analysis

meta <- metagene_profiles(data = serp,
profilefun = make_average_profilefun(data = serp),
len = 500,
bin = “byaa”,
binmethod = “sum”,
align = “start”,
bpparam = BiocParallel::SnowParam(workers = 12))

check

head(meta,3)

# A tibble: 3 × 7

exp rep type id pos summary boot

1 HEK 1 ip 0 0 0.0953 FALSE

2 HEK 1 ip 0 1 0.0867 FALSE

3 HEK 1 ip 0 2 0.0644 FALSE

4×6

plot_metagene_profiles(df = meta,
group = type,
color = type,
# ytrans = “log2”,
ylab = “Ribosome density (AU)”) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
panel.grid = element_blank()) +
scale_fill_manual(values = c(tt = “grey”,ip = “#336699”),
labels = c(“tt” = “Monosome”,”ip” = “Disome”),
name = “”) +
scale_color_manual(values = c(tt = “grey”,ip = “#336699”),
labels = c(“tt” = “Monosome”,”ip” = “Disome”),
name = “”) +
ylim(0,0.3)

❝
使用嘌呤霉素或者蛋白酶抑制剂都会消化共组装形成的新生肽复合物，作者使用着两种药物处理然后进行测序，观察共组装的结果：

我们构建一个新的 serp object：

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

enrichment profile

make serp object for human

serp_new <- load_serp(noPuro =list(disome = c(“GSE151959_RAW/GSM4594351_Dis_HEK_no_puro_Rep1.h5”,
“GSE151959_RAW/GSM4594355_Dis_HEK_no_puro_Rep2.h5”),
mono = c(“GSE151959_RAW/GSM4594350_Mono_HEK_no_puro_Rep1.h5”,
“GSE151959_RAW/GSM4594354_Mono_HEK_no_puro_Rep2.h5”)),
Puro = list(disome = c(“GSE151959_RAW/GSM4594353_Dis_HEK_puro_Rep1.h5”,
“GSE151959_RAW/GSM4594357_Dis_HEK_puro_Rep2.h5”),
mono = c(“GSE151959_RAW/GSM4594352_Mono_HEK_puro_Rep1.h5”,
“GSE151959_RAW/GSM4594356_Mono_HEK_puro_Rep2.h5”)),
no_PK = list(disome = c(“GSE151959_RAW/GSM4594359_Dis_HEK_no_PK.h5”),
mono = c(“GSE151959_RAW/GSM4594358_Mono_HEK_no_PK.h5”)),
low_PK = list(disome = c(“GSE151959_RAW/GSM4594361_Dis_HEK_PK_low.h5”),
mono = c(“GSE151959_RAW/GSM4594360_Mono_HEK_PK_low.h5”)),
mid_PK = list(disome = c(“GSE151959_RAW/GSM4594363_Dis_HEK_PK_mid.h5”),
mono = c(“GSE151959_RAW/GSM4594362_Mono_HEK_PK_mid.h5”)),
high_PK = list(disome = c(“GSE151959_RAW/GSM4594365_Dis_HEK_PK_high.h5”),
mono = c(“GSE151959_RAW/GSM4594364_Mono_HEK_PK_high.h5”)),
veryhigh_PK = list(disome = c(“GSE151959_RAW/GSM4594367_Dis_HEK_PK_veryhigh.h5”),
mono = c(“GSE151959_RAW/GSM4594366_Mono_HEK_PK_veryhigh.h5”)),
bin = ‘byaa’)

save

save(serp_new,file = “serp_new.rda”)
然后分析绘图（怎么作者自己提供的数据画出来还不一样）：

enrichment metagene analysis

meta <- metagene_profiles(data = serp_new,
profilefun = make_enrichment_profilefun(data = serp_new,
sample1 = c(“disome”),
sample2 = c(“mono”)),
len = 500,
bin = “byaa”,
binmethod = “sum”,
align = “start”,
bpparam = BiocParallel::SnowParam(workers = 12))

save

save(meta,file = “meta.rda”)

plot

mete_pk <- plot_metagene_profiles(df = meta %>%
dplyr::filter(exp %in% c(“no_PK”,”low_PK”,”mid_PK”,”high_PK”)),
group = exp,
ylab = “Enrichment (disome / monosome)”,
color = exp) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
panel.grid = element_blank()) +
xlab(“Distance from start (codons)”) +
scale_fill_manual(values = c(no_PK = “black”,low_PK = “#FF9900”,
mid_PK = “#993399”,high_PK = “red”),
name = “”) +
scale_color_manual(values = c(no_PK = “black”,low_PK = “#FF9900”,
mid_PK = “#993399”,high_PK = “red”),
name = “”)

mete_puro <- plot_metagene_profiles(df = meta %>%
dplyr::filter(exp %in% c(“noPuro”,”Puro”)),
group = rep,
ylab = “Enrichment (disome / monosome)”,
color = exp) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
panel.grid = element_blank()) +
xlab(“Distance from start (codons)”) +
scale_fill_manual(values = c(noPuro = “black”,Puro = “#CC3333”),
name = “”) +
scale_color_manual(values = c(noPuro = “black”,Puro = “#CC3333”),
name = “”)

combine 6×6

mete_pk / mete_puro

single gene profile

❝
前面绘制的是 metagene 的，我们画一下文章里的单个基因的 profile：

make serp object for human

serp_human <- load_serp(HEK =list(HEK_disome = c(“GSE151959_RAW/GSM4594341_Dis_HEK_Rep1.h5”,
“GSE151959_RAW/GSM4594343_Dis_HEK_Rep2.h5”),
HEK_mono = c(“GSE151959_RAW/GSM4594340_Mono_HEK_Rep1.h5”,
“GSE151959_RAW/GSM4594342_Mono_HEK_Rep2.h5”),
HEK_no_puro_disome = c(“GSE151959_RAW/GSM4594351_Dis_HEK_no_puro_Rep1.h5”,
“GSE151959_RAW/GSM4594355_Dis_HEK_no_puro_Rep2.h5”),
HEK_no_puro_mono = c(“GSE151959_RAW/GSM4594350_Mono_HEK_no_puro_Rep1.h5”,
“GSE151959_RAW/GSM4594354_Mono_HEK_no_puro_Rep2.h5”),
HEK_puro_disome = c(“GSE151959_RAW/GSM4594353_Dis_HEK_puro_Rep1.h5”,
“GSE151959_RAW/GSM4594357_Dis_HEK_puro_Rep2.h5”),
HEK_puro_mono = c(“GSE151959_RAW/GSM4594352_Mono_HEK_puro_Rep1.h5”,
“GSE151959_RAW/GSM4594356_Mono_HEK_puro_Rep2.h5”),
HEK_no_PK_disome = c(“GSE151959_RAW/GSM4594359_Dis_HEK_no_PK.h5”),
HEK_no_PK_mono = c(“GSE151959_RAW/GSM4594358_Mono_HEK_no_PK.h5”),
HEK_PK_low_disome = c(“GSE151959_RAW/GSM4594361_Dis_HEK_PK_low.h5”),
HEK_PK_low_mono = c(“GSE151959_RAW/GSM4594360_Mono_HEK_PK_low.h5”),
HEK_PK_mid_disome = c(“GSE151959_RAW/GSM4594363_Dis_HEK_PK_mid.h5”),
HEK_PK_mid_mono = c(“GSE151959_RAW/GSM4594362_Mono_HEK_PK_mid.h5”),
HEK_PK_high_disome = c(“GSE151959_RAW/GSM4594365_Dis_HEK_PK_high.h5”),
HEK_PK_high_mono = c(“GSE151959_RAW/GSM4594364_Mono_HEK_PK_high.h5”),
HEK_PK_veryhigh_disome = c(“GSE151959_RAW/GSM4594367_Dis_HEK_PK_veryhigh.h5”),
HEK_PK_veryhigh_mono = c(“GSE151959_RAW/GSM4594366_Mono_HEK_PK_veryhigh.h5”)),
U2OS = list(U2OS_disome = c(“GSE151959_RAW/GSM4594345_Dis_U2OS_Rep1.h5”,
“GSE151959_RAW/GSM4594347_Dis_U2OS_Rep2.h5”),
U2OS_mono = c(“GSE151959_RAW/GSM4594344_Mono_U2OS_Rep1.h5”,
“GSE151959_RAW/GSM4594346_Mono_U2OS_Rep2.h5”)),
bin = ‘byaa’)

save(serp_human,file = “serp_human.rda”)

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

plot for each gene

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

DCTN1_hek <-
plot(data = serp_human,
gene = “DCTN1”,
type = “rpm”,
samples = c(“HEK_mono”,”HEK_disome”),
exp = “HEK”,
bin = “byaa”,
window_size = 45,
conf.level = 0.95,
color = sample) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
panel.grid = element_blank()) +
scale_fill_manual(values = c(HEK_mono = “black”,HEK_disome = “#003366”),
name = “”) +
guides(alpha = F) +
xlab(“Position in CDS (codons)”) +
ylab(“Ribosome density (RPM)”) +
scale_y_continuous()

DCTN1_u2os <-
plot(data = serp_human,
gene = “DCTN1”,
type = “rpm”,
samples = c(“U2OS_mono”,”U2OS_disome”),
exp = “U2OS”,
bin = “byaa”,
window_size = 45,
conf.level = 0.95,
color = sample) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
panel.grid = element_blank()) +
scale_fill_manual(values = c(U2OS_mono = “black”,U2OS_disome = “#003366”),
name = “”) +
guides(alpha = F) +
xlab(“Position in CDS (codons)”) +
ylab(“Ribosome density (RPM)”) +
scale_y_continuous()

library(patchwork)

5×5

DCTN1_hek/DCTN1_u2os

❝
前面两个基因是没有共组装形成的，后面两个基因则是存在共组装现象的：

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

g1 <- c(“JUN”,”SRP54″,”CAPRIN1″,”NFKB1″)

loop plot

x = 1

lapply(seq_along(g1),function(x){
phek <-
plot(data = serp_human,
gene = g1[x],
type = “rpm”,
samples = c(“HEK_mono”,”HEK_disome”),
exp = “HEK”,
bin = “byaa”,
window_size = 45,
conf.level = 0.95,
color = sample) +
labs(title = paste(g1[x],”(HEK-293 T cells)”)) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
plot.title = element_text(hjust = 0.5,face = “bold.italic”,size = rel(1)),
panel.grid = element_blank()) +
scale_fill_manual(values = c(HEK_mono = “black”,HEK_disome = “#003366”),
name = “”) +
guides(alpha = F,fill = F) +
xlab(“Position in CDS (codons)”) +
ylab(“Ribosome density (RPM)”) +
scale_y_continuous() + xlab(“”)

pu2os <-
plot(data = serp_human,
gene = g1[x],
type = “rpm”,
samples = c(“U2OS_mono”,”U2OS_disome”),
exp = “U2OS”,
bin = “byaa”,
window_size = 45,
conf.level = 0.95,
color = sample) +
labs(title = paste(g1[x],”(U2OS cells)”)) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
plot.title = element_text(hjust = 0.5,face = “bold.italic”,size = rel(1)),
panel.grid = element_blank()) +
scale_fill_manual(values = c(U2OS_mono = “black”,U2OS_disome = “#003366”),
name = “”) +
guides(alpha = F,fill = F) +
xlab(“Position in CDS (codons)”) +
ylab(“Ribosome density (RPM)”) +
scale_y_continuous()

cowplot::plot_grid(plotlist = list(phek,pu2os),ncol = 1)
}) -> plist

combine plots 6×15

cowplot::plot_grid(plotlist = plist,nrow = 1)

❝
不同药物处理对共组装结果的影响：

plot for pk

pk <-
plot(data = serp_new,
gene = “DCTN1”,
type = “enrichment”,
sample1 = c(“disome”),
sample2 = c(“mono”),
exp = c(“no_PK”,”low_PK”,”mid_PK”,”high_PK”),
bin = “byaa”,
window_size = 45,
conf.level = 0.95) +
labs(title = paste(“DCTN1″,”(HEK-293 T cells)”)) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
plot.title = element_text(hjust = 0.5,face = “bold.italic”,size = rel(1)),
panel.grid = element_blank()) +
scale_fill_manual(values = c(no_PK = “black”,low_PK = “#FF9900”,
mid_PK = “#993399”,high_PK = “red”),
name = “”) +
scale_color_manual(values = c(no_PK = “black”,low_PK = “#FF9900”,
mid_PK = “#993399”,high_PK = “red”),
name = “”) +
guides(alpha = F) +
xlab(“Position in CDS (codons)”) +
ylab(“Enrichment (disome / monosome)”)

plot for puro

puro <-
plot(data = serp_new,
gene = “DCTN1”,
type = “enrichment”,
sample1 = c(“disome”),
sample2 = c(“mono”),
exp = c(“noPuro”,”Puro”),
bin = “byaa”,
window_size = 45,
conf.level = 0.95) +
labs(title = paste(“DCTN1″,”(HEK-293 T cells)”)) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
plot.title = element_text(hjust = 0.5,face = “bold.italic”,size = rel(1)),
panel.grid = element_blank()) +
scale_fill_manual(values = c(noPuro = “black”,Puro = “#CC3333”),
name = “”) +
scale_color_manual(values = c(noPuro = “black”,Puro = “#CC3333”),
name = “”) +
guides(alpha = F) +
xlab(“Position in CDS (codons)”) +
ylab(“Enrichment (disome / monosome)”)

6×6

pk/puro

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

g1 <- c(“JUN”,”SRP54″,”CAPRIN1″,”NFKB1″)

loop plot

x = 1
lapply(seq_along(g1),function(x){
pk <-
plot(data = serp_new,
gene = g1[x],
type = “enrichment”,
sample1 = c(“disome”),
sample2 = c(“mono”),
exp = c(“no_PK”,”low_PK”,”mid_PK”,”high_PK”),
bin = “byaa”,
window_size = 45,
conf.level = 0.95) +
labs(title = paste(g1[x],”(HEK-293 T cells)”)) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
legend.position = “none”,
plot.title = element_text(hjust = 0.5,face = “bold.italic”,size = rel(1)),
panel.grid = element_blank()) +
scale_fill_manual(values = c(no_PK = “black”,low_PK = “#FF9900”,
mid_PK = “#993399”,high_PK = “red”),
name = “”) +
scale_color_manual(values = c(no_PK = “black”,low_PK = “#FF9900”,
mid_PK = “#993399”,high_PK = “red”),
name = “”) +
guides(alpha = F) +
# xlab(“Position in CDS (codons)”) +
xlab(“”) +
ylab(“Enrichment (disome / monosome)”)

# plot for puro
puro <-
plot(data = serp_new,
gene = g1[x],
type = “enrichment”,
sample1 = c(“disome”),
sample2 = c(“mono”),
exp = c(“noPuro”,”Puro”),
bin = “byaa”,
window_size = 45,
conf.level = 0.95) +
theme_bw() +
theme(axis.text = element_text(colour = “black”),
legend.position = “none”,
plot.title = element_blank(),
panel.grid = element_blank()) +
scale_fill_manual(values = c(noPuro = “black”,Puro = “#CC3333”),
name = “”) +
scale_color_manual(values = c(noPuro = “black”,Puro = “#CC3333”),
name = “”) +
guides(alpha = F) +
xlab(“Position in CDS (codons)”) +
ylab(“Enrichment (disome / monosome)”)

return(cowplot::plot_grid(plotlist = list(pk,puro),ncol = 1))
}) -> plist2

combine plots 6×15

cowplot::plot_grid(plotlist = plist2,nrow = 1)

双曲线模型选择高置信度共组装靶标

❝
存在翻译共组装的基因的 enrichment（disome/monosome） 的 profile 应该是 S 型曲线，组装提出了两种共组装翻译的模式(顺式和反式)，它们所对应的组装的情况和曲线也相应不同。为了筛选高置信度的候选基因，作者规定了 5 个规则 来进行筛选。

S 型曲线建模（3 个模型）

定义高置信度候选基因

总结

❝
只是复现了一小部分东西，感觉文章的很多东西都超越了自己的水平，很多东西实现不了，觉得自己学识还是太浅薄了，只能望洋兴叹了，哈哈。
结尾

❝
路漫漫其修远兮,吾将上下而求索。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。