中文
PK（藥代動力學）分析主要研究藥物在體內的吸收、分布、代謝和排泄過程，其涉及的相關圖表包括血藥濃度-時間曲線、PK參數計算、劑量線性分析等。
R 是一款開源免費、擴展性強、腳本化操作的統計編程語言和環境，專為數據分析和可視化設計，適用于早期研究、探索性分析、內部決策支持等非遞交階段。正式遞交仍需使用已驗證的商業軟件（如Phoenix WinNonlin、SAS），以確保合規性。
本文演示使用R完成PK分析全流程，包含：
?
非房室模型（NCA）：用于計算關鍵PK參數（如Cmax、AUC），無需假設房室結構，適用于早期藥物研發階段。
?
Power Model分析：探索劑量與暴露量的線性關系，支持劑量優化決策。
?
可視化分析：通過專業圖表（如線性和半對數藥時曲線圖）直觀呈現藥物動力學特征。
?
自動化報告生成：利用r2rtf包實現從數據到Word報告的一鍵輸出，減少人工操作誤差。
流程結合了開源工具（如PKNCA、ggplot2）與傳統商業軟件（WinNonlin、SAS）的對比驗證，確保分析結果的可靠性。
# 核心分析包library(tidyverse)library(readxl)library(PKNCA) # 非房室分析library(r2rtf) # RTF輸出# 輔助工具library(knitr) # 交互式表格
PART 01
數據模擬
原始數據模擬
pc <- read_xlsx(path=str_c(getwd(),"/outputs/pc.xlsx"))# 查看模擬數據knitr::kable(pc[1:10, 1:6])
PART 02
非房室分析（NCA）
關鍵PK參數計算
# 設置PK參數計算的單位轉換表ppunits <- pknca_units_table(concu="ug/mL", doseu="mg/kg", amountu="mg", timeu="hr")conc_raw <- arrange(pc, SUBJID,ATPTN)# 設置濃度數據conc_obj_multi <- PKNCAconc(data=conc_raw,formula=AVAL~ATPTN|SUBJID)# 設置劑量數據dose_raw <- select(conc_raw,ATPTN,DOSE,SUBJID) %>% distinct()dose_obj_multi <- PKNCAdose(dose_raw,DOSE~ATPTN|SUBJID,route="intravascular",duration=1.5)# 設置參數計算時間區間intervals_manual <- data.frame(start=0,end=504,cmax=TRUE,auclast=TRUE,aucinf.obs=TRUE)o_data <- PKNCAdata(conc_obj_multi, dose_obj_multi,intervals=intervals_manual,units=ppunits)# 計算PK參數o_results <- pk.nca(o_data)# 結果呈現ncasingle <- o_results[["result"]] %>%select(SUBJID,start,end,PPTESTCD,PPORRES,PPORRESU) %>%filter(PPTESTCD %in% c("cmax","auclast","aucinf.obs"))# 輸出PK參數kable(pivot_wider(ncasingle,id_cols=1,names_from = "PPTESTCD", values_from = "PPORRES"))
WinNonlin計算結果
PART 03
Power Model分析
劑量-暴露量關系建模
Power Model (原始形式):
線性化形式 (取自然對數后):
·
數據轉換：對劑量和PK參數取自然對數，轉化為線性模型。
·
擬合與檢驗：使用lm函數進行回歸分析，計算斜率置信區間和R2，評估模型擬合優度。
· 可視化：展示不同參數的劑量效應趨勢，圖中標注斜率、R2及90%置信區間，直觀支持統計結論。
# 劑量與參數的自然對數rst <- left_join(ncasingle,distinct(dose_raw,DOSE,SUBJID),by="SUBJID") %>%mutate(logx=log(DOSE),logy=log(PPORRES),PPTESTCD=factor(PPTESTCD,levels=c("cmax","auclast","aucinf.obs")))# 模型擬合split <- split(rst,~PPTESTCD,drop=TRUE)lmci <- map(split,\(x) data.frame(PPTESTCD=unique(x$PPTESTCD),logx=lm(logy~logx,x)[["coefficients"]][["logx"]],Rsq=summary(lm(logy~logx,x))[["r.squared"]],confint(lm(logy~logx,x), "logx", level= 0.90 ))) %>%bind_rows() %>%# 設置繪圖呈現格式mutate(.logx=formatC(`logx`,digits=4,format = "f"),lwci=formatC(`X5..`,digits=4,format = "f"),upci=formatC(`X95..`,digits=4,format = "f"),Rsq=formatC(`Rsq`,digits=4,format = "f"),ci=str_c("Rsq:",Rsq,";Slope:",.logx,";\n90%CI (",lwci,", ",upci,")"),) %>%select(PPTESTCD,ci)# 劑量線性繪圖呈現(base <- ggplot(rst,aes(x=logx,y=logy)) +geom_point(shape=24) +geom_smooth(method = "lm", se = TRUE) +geom_text(aes(x=1.5,y=-Inf,label=ci),data=lmci,size=3,vjust=-0.5) +scale_x_continuous(breaks=c(0.69,1.38,2.08),labels=c("Ln(2)","Ln(4)","Ln(8)")) +labs(title = "Power Model",x="ln(Dose)(mg/kg)",y="ln(param)") +facet_wrap(vars(PPTESTCD),scales = "free", axes="all" ,labeller = as_labeller(c("cmax"="C[max](ug/mL)","auclast"="AUC[last](h%.%ug/mL)","aucinf.obs"="AUC[0-inf](h%.%ug/mL)"),,default = label_parsed)) +theme_minimal())
SAS 結果
通過SAS代碼復現分析，結果一致性驗證了R流程的可靠性。
PART 04
可視化分析
平均血藥濃度曲線
·
線性尺度：展示各劑量組濃度隨時間的變化趨勢，誤差條表示標準差，反映個體間變異。
· 對數尺度：通過scale_y_log10轉換，突出低濃度階段的細節，適用于寬動態范圍的數據。
pcsas % mutate( aval=as.numeric(AVAL), DOSE=str_c(DOSE,DOSEU) ) # 分組時點血藥濃度定量統計stplot % group_by(DOSE,ATPTN) %>% summarise( Mean=mean(aval,na.rm = TRUE), SD=sd(aval,na.rm = TRUE), lower=ifelse(is.na(SD) | Mean-SD % ungroup()# 線性血藥濃度-時間曲線(normal % ggplot(aes(x = ATPTN, y = Mean ,colour=DOSE)) + geom_errorbar(mapping=aes(ymin=lower, ymax=upper),width = 0.5,alpha=0.5) + geom_line(linewidth = 0.5,alpha=0.5) + geom_point(alpha=0.5) + scale_x_continuous(breaks=c(24,72,168,336,504)) + labs(x="Time(h)",y="Concentration(mean)",colour="劑量組:",title = "藥時曲線") + theme_classic() )
圖 藥時曲線
# 半對數血藥濃度-時間曲線normal %+% subset(stplot, Mean>0) + scale_y_log10(n.breaks=4,labels = scales::trans_format("log10", scales::math_format(10^.x))) + annotation_logticks(sides = "l",outside = TRUE) + expand_limits(y = c(1,10,100)) + coord_cartesian(clip = "off")
圖 藥時曲線
PART 05
報告生成
Word文檔自動導出
工具鏈：r2rtf包支持從數據匯總到RTF格式報告的全流程自動化，功能類似SAS的Proc Report。
關鍵步驟：
·
數據匯總：按劑量組和時間點計算均值、標準差、幾何均值等統計量。
·
表格格式化：通過pivot_wider和arsenal包生成結構化表格。
·
RTF輸出：定制頁面方向（橫向）、邊距和分欄布局，生成格式精美、便于審閱的標準化報告。
pcsas <- pc %>%mutate(aval=as.numeric(AVAL),log=ifelse(aval>0,log(aval),NA),order=1,DOSE=str_c(DOSE,DOSEU),PCSTRESC=str_trim(formatC(AVAL,digits = 3,format = "f")) ,)# 定量統計，并調整樣式st <- filter(pcsas,!is.na(aval)) %>%group_by(DOSE,ATPTN) %>%summarise(N=paste(sum(!is.na(aval))),Mean=str_trim(formatC(mean(aval,na.rm = TRUE),digits = 3,format = "f")),SD=str_trim(formatC(sd(aval,na.rm = TRUE),digits = 3,format = "f")) ,Min=formatC(min(aval,na.rm = TRUE),digits = 2,format = "f"),Median=str_trim(formatC(median(aval,na.rm = TRUE),digits = 2,format = "f")),Max=formatC(max(aval,na.rm = TRUE),digits = 2,format = "f"),`CV%`=case_when(mean(aval,na.rm = TRUE) > 0 ~ formatC(sd(aval,na.rm = TRUE)/mean(aval,na.rm = TRUE)*100,digits = 1,format = "f"),.default = NA),Geomean=case_when(mean(aval,na.rm = TRUE) > 0 ~ str_trim(formatC(exp(mean(log,na.rm = TRUE)),digits = 3,format = "f")),.default = NA),) %>%pivot_longer(cols = c("N","Mean","SD","Min","Median","Max","CV%","Geomean"),names_to = "SUBJID",values_to = "PCSTRESC") %>%mutate(order=case_when(SUBJID=="N" ~ 2 ,SUBJID=="Mean" ~ 3 ,SUBJID=="SD" ~ 4 ,SUBJID=="Min" ~ 5 ,SUBJID=="Median" ~ 6 ,SUBJID=="Max" ~ 7 ,SUBJID=="CV%" ~ 8 ,SUBJID=="Geomean" ~ 9),SUBJID=str_c("\\b1{",SUBJID,"}"),PCSTRESC=str_c("\\b{",PCSTRESC,"}"))# 合并原始數據與統計結果rst <- bind_rows(pcsas,st) %>%select(DOSE,order,SUBJID,ATPTN,PCSTRESC) %>%arrange(DOSE,order,SUBJID,ATPTN) %>%pivot_wider(id_cols = c(1:3),names_from = c("ATPTN"),values_from = "PCSTRESC") %>%ungroup() %>%select(-order) %>%group_by(DOSE) %>%group_modify(~ add_row(.x))# 設置輸出樣式rtf <- rst %>%rtf_page(orientation = "landscape",border_first = "single",border_last = "single",margin = c(0.75, 0.75, 0.75, 0.75, 0.5, 0.5),nrow=100) %>%rtf_colheader("||Time",border_left=NULL,border_right = NULL,border_top = "single",border_bottom = c("","","single"),col_rel_width = c(8,8,96)) %>%rtf_colheader(colheader = "劑量|受試者|Pre|EOI|0.5h|1h|2h|6h|D2|D3|D4|D8|D15|D22-Pre",border_left=NULL,border_right = NULL,border_top = NULL,border_bottom = "single",,col_rel_width = rep(9,14)) %>%rtf_body(border_left = NULL,border_right = NULL,pageby_row = "column",group_by="DOSE",subline_by=NULL) %>%rtf_encode() %>%str_replace('fcharset161','fcharset134') %>%write_rtf(str_c("outputs/","table_pc_summary.rtf"))
表 模擬數據結果
PART 06
總結
本文通過R語言實現了藥代動力學的全流程分析，涵蓋數據預處理、參數計算、統計建模、可視化及報告生成，適用于非遞交的PK數據分析。熙寧生物憑借專業的PK編程團隊和資深藥理專家，結合驗證工具（如Winnonlin、SAS）與開源軟件（如R），能夠高效精準地滿足客戶各類PK/PD統計分析需求，為藥物研發與臨床研究提供可靠支持。
PART 07
熙寧生物臨床藥理服務平臺
熙寧生物臨床藥理服務平臺，專注于創新藥（涵蓋生物大分子及小分子）的臨床藥理學研究，擁有豐富的項目經驗和專業的技術能力。我們提供全方位的臨床藥理學服務，包括：
?
PK NCA分析：采用WinNonlin軟件非房室模型（NCA）計算PK參數，結合SAS和R等專業軟件進行數據編程，確保分析結果精準可靠。
?
實時分析支持劑量爬坡：通過實時數據分析，為SMC會議提供科學依據，優化臨床試驗設計。
?
臨床藥理學研究設計與方案撰寫：基于客戶需求，定制化設計研究方案，確?？茖W性和合規性。
?
完整PK/PD統計分析及CSR撰寫：提供生物分析檢測到統計分析報告生成的一站式服務，確保高效交付。
平臺配備資深的臨床藥理專家、專業統計師及統計編程團隊，能夠對統計分析報告中的數據進行深度解讀，提供專業的洞見和建議。我們嚴格按照CDISC標準進行統計編程，確保交付成果符合NMPA、FDA等國際監管機構的要求，助力藥物快速通過審批。
PART 08
我們的優勢
● 在生物大分子領域，我們擁有豐富的PK統計分析經驗，涵蓋 CAR-T細胞療法、單抗、雙抗、ADC（抗體藥物偶聯物）等藥物類型，能夠精準分析其非線性動力學特征、靶點介導的藥物處置（TMDD）等特殊機制。
● 對小分子藥物，我們深入分析其吸收、分布、代謝和排泄（ADME）特性，全面評估其藥代動力學行為。
熟練使用 SAS 和 R 進行數據編程，生成高度定制化的高質量圖表，相比 WinNonlin 的基礎功能，R 和 SAS 在靈活性和可視化效果上更具優勢，確保數據呈現清晰直觀。
憑借高效的團隊協作和先進的技術工具，我們能夠快速響應客戶需求，提供高效、精準的交付成果，助力藥物研發加速推進。
參考文獻：
[1] The PKNCA R Package
[2] Introduction to r2rtf
[3] Phoenix WinNonlin User’s Guide
