越來越懶了隔缀，不想寫記錄呀幔妨！算了鹦赎，強迫自己一把！
水平有限误堡，若有錯誤古话，懇請指出！
信使RNA ( messenger RNA, mRNA)加工過程中的精細調(diào)控對基因的表達具有重要影響,是產(chǎn)生基因功能多樣性的重要機制锁施。在前體mRNA成熟過程中,環(huán)境的細微變化能夠?qū)е略?b>mRNA的不同剪接位點上進行選擇性的剪接和多聚腺苷酸化,此過程被稱為可選擇性多聚腺苷酸化(alerna"tive polyadenylation , APA)陪踩。可變多聚腺苷酸化(alternative polyadenylation, APA)是真核細胞mRNA成熟過程中針對前體mRNA 3′端的一種加工修飾方式,是重要的轉(zhuǎn)錄后調(diào)控機制悉抵。
近年來,APA作為調(diào)控基因表達的重要方式越來越受到重視肩狂。APA現(xiàn)象廣泛存在于真核生物中,絕大多數(shù)真核基因具有多個多聚腺苷酸[poly( A)]位點。例如,人類有70% ~75%的基因存在APA事件,果蠅1姥饰、蠕蟲和斑馬魚[3也有約半數(shù)的基因存在這類現(xiàn)象傻谁。APA能產(chǎn)生長短不同的mRNA轉(zhuǎn)錄本,不同的mRNA轉(zhuǎn)錄本具有不同的穩(wěn)定性、翻譯能力和表達水平列粪，從而產(chǎn)生多樣化的生物學功能审磁。
?APA的形成過程APA形成過程包括poly (A)位點的選擇和多聚腺苷酸化,這兩個過程十分重要并受到多種因素的調(diào)控。選擇不同的多聚腺苷酸化信號是形成APA的關(guān)鍵篱竭。
經(jīng)典的poly (A)模體是AAUAAA,其他非經(jīng)典的poly (A)模體變體包括. AUUAAA, AAGAAA和UAUAAA等力图。在人類基因組中,AAUAAA約占所有多聚腺苷酸化模體的50%,AUUAAA約占16%、A[A/U]UAAA中約20%存在單核苷酸變異,只有約10%的多聚腺苷酸化模體沒有類似于AAUAAA的識別區(qū)域5mRNA的多聚腺苷酸化包含2個必要步驟:1前體mRNA在poly(A)位點進行剪接;2在3'端加poly (A)尾掺逼。
多聚腺苷酸化的主要功能是保證轉(zhuǎn)錄本的高效輸出吃媒、翻譯以及穩(wěn)定。多聚腺苷酸化的順利形成是前體mRNA序列中的順式作用元件、反式作用因子以及細胞核中多種酶和蛋白因子相互作用的結(jié)果赘那。
選擇性多聚腺苷酸化（APA）是一種重要的前體RNA加工機制刑桑，廣泛存在于所有真核生物中。通過在RNA 3′UTR不同位置上添加polyA尾巴募舟，可以選擇性地調(diào)節(jié)3′UTR的長短祠斧。由于3′UTR含有多種順式調(diào)控元件，例如：miRNA或RNA結(jié)合蛋白（RBP）結(jié)合位點拱礁，因此琢锋，APA可以通過調(diào)控?3′非翻譯區(qū)(3′UTR)?的長度，影響目標mRNA的穩(wěn)定性和翻譯效率以及翻譯后蛋白質(zhì)的細胞定位呢灶，進而精細調(diào)節(jié)基因表達吴超，對一系列細胞過程（如增殖、分化和腫瘤發(fā)生）產(chǎn)生根本性的影響鸯乃。
內(nèi)含子多聚腺苷酸化(intronic polyadenylation, IPA)通過形成丟失重要結(jié)構(gòu)域的截短型蛋白實現(xiàn)對靶基因的調(diào)控,參與形成腫瘤新抗原。APA具有腫瘤特異性,有可能用于腫瘤分子分型和靶向治療缨睡。
可變多聚腺苷酸化Alternative?Polyadenylation (APA)鸟悴，如下圖所示（圖片來自參考），在不同的APA信號位點切割奖年，然后添加polyA细诸。這種調(diào)控機制屬于轉(zhuǎn)錄后調(diào)控，可能會影響蛋白的序列（發(fā)生在編碼區(qū)）拾并，也可能影響蛋白的穩(wěn)定性（比如非編碼區(qū)內(nèi)的miRNA的調(diào)控區(qū)域）揍堰。其實也是可變剪接的一種情況。
一般APA可以分為四種類型：
1 嗅义、3’UTR APA：發(fā)生在末端外顯子內(nèi)屏歹，產(chǎn)生具有不同長度3’UTR的轉(zhuǎn)錄本，不影響蛋白編碼功能之碗，是最常見的APA形式蝙眶；
2、 可變末端外顯子APA：這種APA產(chǎn)生了末端外顯子不同的轉(zhuǎn)錄本褪那，影響蛋白編碼功能幽纷；
3 、內(nèi)含子APA：發(fā)生于在內(nèi)含子區(qū)博敬，延長了某個內(nèi)部外顯子并使之成為末端外顯子友浸；
4、 內(nèi)部外顯子APA：在編碼區(qū)域內(nèi)部發(fā)生剪切和多聚腺苷酸化偏窝。
常用的軟件是Dapars收恢，這個軟件現(xiàn)在也有了升級的版本Dapars2武学。參考：
https://github.com/ZhengXia/dapars
https://github.com/3UTR/DaPars2
分析流程很相似，Dapars2多了 normalize library sizes 伦意。
安裝軟件
官網(wǎng)：https://github.com/3UTR/DaPars2
source activate dapars2? ?# 需要python2.7的環(huán)境火窒！
cd? /data/APA_DaPars2? ?# 分析路徑
step1 :這一步有兩個文件需要從UCSC下載：whole_gene.bed,? id_from_UCSC.txt。
1驮肉、下載：hg38_refseq_whole_gene.bed
genome: Human
assembly:hg38
group: Genes and Gene Predictions
track: NCBI_REfSeq
table: refGene All
region: genome
output format: BED - browser extensible data
output file: hg38_refseq_whole_gene.bed
點‘get output’ button,下一頁點‘Output refGene as BED’ 再點 ‘get output’ button.
2熏矿、下載：*hg38_Refseq_id_from_UCSC.txt
genome: Human
assembly: hg38
group: Genes and Gene Predictions
track: NCBI REfSeq
table: refGene? All
region: genome
output format: selected fields from primary and related tables
output file:?hg38_Refseq_id_from_UCSC.txt
點 ‘get output’ button，下一個界面選擇：
name: Name of gene (usually transcript_id from GTF)
name2: Alternate name (e.g. gene_id from GTF)
點 ‘get output’ 保存文件
head? hg38_refseq_whole_gene.bed? ? ? hg38_Refseq_id_from_UCSC.tx
DaPars_Extract_Anno.py? 腳本 構(gòu)建注釋
python DaPars_Extract_Anno.py -b hg38_refseq_whole_gene.bed? -s? hg38_Refseq_id_from_UCSC.txt? -o hg38_refseq_extracted_3UTR.bed
##### step2 生成 bedgraph
bedtools genomecov -ibam? test.sorted.bam? -bga -split -trackline? >? test.sorted.bam.out.wig
### 去掉wig的第一行离钝，并加上chr?
ls *wig | while read a;do sed -r -i -e '1d' -e 's/^(.?.\t)/chr\1/' $a;done &
head??test.sorted.bam.out.wig
### 看下共有哪些染色體及覆蓋度
cut -f1 test.sorted.bam.out.wig? | uniq? -c
挑選出需要的chr染色體
ls *wig | while read a ; do grep "chr" $a > ../wig/$a ; done &
?step3 生成 mapping_bam_location_with_depth.txt, 拿來做測序深度歸一化
## 我用 samtools? ?flagstat 統(tǒng)計bam 文件 的 mapping reads票编，all.samtools_flagstat_result.xls
cat all.samtools_flagstat_result.xls
## cut? 出? ? mapping reads，?第四部的?configure文件需要,???需和wig文件名奈辰、順序一致栏妖。
cut -f1,6 all.samtools_flagstat_result.xls |sed 's#(.*##g' | sed? 's#R#R.sorted.bam.out.wig#'? > all.samtools_flagstat.mapped.txt
head? ? ?all.samtools_flagstat.mapped.txt
step4? 寫configure_file
###? 列出wig文件奖恰，以逗號分割。批量操作宛裕，不想手打瑟啃， 粘貼在下面的 configure文件里
cut -f1 all.samtools_flagstat.mapped.txt | sed 's#^#./wig/#g'| paste -s -d ","
# 寫configure_file
cat DaPars2_test_data_configure.txt
Annotated_3UTR=?hg38_refseq_extracted_3UTR.bed
Aligned_Wig_files= sample1.wig, sample2.wig #以逗號分割 ，列出wig文件
Output_directory= result_DaPars2_data/
Output_result_file= result_DaPars2_data
sequencing_depth_file=?all.samtools_flagstat.mapped.txt
Coverage_threshold=10
Num_threads=30? ## 30個線程
# 運行? DaPars2.sh
python? ?DaPars2_Multi_Sample_Multi_Chr.py? DaPars2_test_data_configure.txt
###?如果要單獨分析每個染色體以節(jié)省時間和內(nèi)存請運行以下命令揩尸，可以每條染色體寫個腳本蛹屿，同時提交，同時運行岩榆，節(jié)省時間和內(nèi)存 错负！
python? ? Dapars2_Multi_Sample.py? DaPars2_test_data_configure.txt? ? chr3
## 看 看chr3
head? result_DaPars2_Test_data_chr3/*txt
Predicted_Proximal_APA' : 表示預測的近端 Poly(A) 位點，'Loci' 表示轉(zhuǎn)錄本的 3'UTR 區(qū)域勇边，其他列代表每個樣本的 PDUI 值犹撒。
再把所有chr 染色體合并到一個文件? all.chr.DaPars2.txt?里。
head??all.chr.DaPars2.txt
我有個問題呀粒褒，如果我想將 實驗組 與 對照組 進行差異分析呢识颊？ 需要怎么實現(xiàn)？我在DaPars2 更新后的版本里奕坟，我沒有看到這種用法呢祥款？在上一步的configure_file里， 我是將 N 和 T 的wig放在一組運行的月杉。 但怎么做腫瘤組和正常組的差異APA分析嗎刃跛？? ?再研究研究~~
后續(xù)來了： 我去給DaPars2作者留言了，作者的回復是：?DaPars2不能做差異分析苛萎，還是建議用以前的DaPars~~~~桨昙， 好吧跌帐，下面重頭開始運行DaPars！
以下是 運行第一版的 DaPars記錄
https://github.com/ZhengXia/dapars
安裝第一版的 DaPars绊率，需 注意 pip install rpy2==2.8.6
第一版的 DaPars版本 的 config 文件：?DaPars_configure
Annotated_3UTR=?hg38_refseq_extracted_3UTR.bed
Group1_aligned_Wig= T1.wig谨敛，T2.wig，T3wig滤否，T4.wig
Group2_aligned_Wig=N1.wig脸狸，?N2.wig，?N3wig藐俺，?N4.wig
Num_least_in_group1=1
Num_least_in_group2=1
Coverage_cutoff=30?
FDR_cutoff=0.05?
PDUI_cutoff=0.4?
Fold_change_cutoff=0.59
運行腳本DaPars_main.py
python? ? ? ?/data/sorftware/dapars-master/src/DaPars_main.py? ? ?DaPars_configure
運行?DaPars_main.py 后的結(jié)果：
在R中炊甲，我們?nèi)〕龅谝涣術(shù)ene，以及后5列的結(jié)果：
library(dplyr)
a = read.table("test.new_DaPars_Test_data_All_Prediction_Results.txt", header=T, row.names= 1, sep="\t", quote="")
dim(a) ; head(a[, (ncol(a)-5): ncol(a)])
b =? ?a[, (ncol(a)-5) :ncol(a)]
b =? ?b[,c(1,2,3,5)];head(b)
b =? b %>% dplyr::rename("TumorMeanPDUI"="Group_A_Mean_PDUI", "NormalMeanPDUI" = "Group_B_Mean_PDUI");head(b)
write.table(b, "PDUI-Pval.txt",sep="\t", quote=F)
system(paste("sed -i '1 s/^/Gene\t/'", "PDUI-Pval.txt"))
system("less -S? PDUI-Pval.txt | sed -e '1s/^Gene\t/gene_id\tgene_name\tchr\tstrand\t/' -e 's/|/\t/g' | cut -f1,2,5,6,7,8? > new.PDUI-Pval.txt")
處理后欲芹，拿來畫圖的new.PDUI-Pval.txt：
R畫圖：
PDUIPVal? <-? read.table("new.PDUI-Pval.txt",row.names=1,header=T)
head(PDUIPVal)
## 1、以0.05位閾值挑出差異APA的基因菱父，
PDUIPVal[which(abs(PDUIPVal[,4]) <=0.05),'DEG0.05'] <- 'no diff'?
PDUIPVal[which((abs(PDUIPVal[,4])>=0.05)&(PDUIPVal[,3]>PDUIPVal[,2])),'DEG0.05'] <- 'Shortened(0.05)'
PDUIPVal[which((abs(PDUIPVal[,4])>=0.05)&(PDUIPVal[,3]<PDUIPVal[,2])),'DEG0.05'] <- 'Lengthened(0.05)'
table(PDUIPVal$DEG0.05)
## 2颈娜、或者要求更嚴格，以0.1位閾值挑出差異APA的基因浙宜，
PDUIPVal[which(abs(PDUIPVal[,4]) <=0.1),'DEG0.1'] <- 'no diff'
PDUIPVal[which((abs(PDUIPVal[,4])>=0.1)&(PDUIPVal[,3]>PDUIPVal[,2])),'DEG0.1'] <- 'Shortened(0.1)'
PDUIPVal[which((abs(PDUIPVal[,4])>=0.1)&(PDUIPVal[,3]<PDUIPVal[,2])),'DEG0.1'] <- 'Lengthened(0.1)'
table(PDUIPVal$DEG0.1)
head(PDUIPVal)
write.table(PDUIPVal, "new.PDUI-Pval.txt", row.names=T, col.names=T, sep="\t",quote=F)
system(paste("sed -i '1 s/^/gene_id\t/'", "new.PDUI-Pval.txt"))
PDUIPVal <- read.table("new.PDUI-Pval.txt",row.names=1,header=T,sep="\t",quote="")
head(PDUIPVal);table(PDUIPVal$DEG0.05);table(PDUIPVal$DEG0.1)
## Figure 1 notes官辽，以0.05為 閾值
pdf("fig1-pduipaad.pdf",width=13,height=5.7)
? par(mfrow=c(1,2))
? par(cex.lab=1.3,cex.axis=1.5)
? par(mar=c(5,5,3,3))
? col0<-rep("#00000022",nrow(PDUIPVal))
? col0[(abs(PDUIPVal[,4])>= 0.05)&(PDUIPVal[,3]>PDUIPVal[,2])] <- "#FF000090" #紅色
? col0[(abs(PDUIPVal[,4])>= 0.05)&(PDUIPVal[,3]<PDUIPVal[,2])]<- "#0000FF90" #藍色
? pch0<-rep(21,nrow(PDUIPVal))
? par(pty="s")
? ? plot(PDUIPVal[,3:2],bg=col0,pch=pch0,col=NA,cex=1.5,axes=FALSE,xlab="",ylab="")
? ? title(xlab="PDUI Normal",ylab="PDUI Tumor",line= 2.3,cex.lab= 1.3)
? ? title(main="",line= 0.5)
? ? box();axis(1,at=c(0,.5,1));axis(2,at=c(0,.5,1))
? ? legend("topleft", c("Lengthened", "Shortened"), pt.bg=c("#0000FF90","#FF000090"), pch=c(21,21), cex=1)
? ? abline(c(.05,1),lty=2,lwd=2)
? ? abline(c(-.05,1),lty=2,lwd=2)
? par(pty="m")
? ? plot((PDUIPVal[,4]),-log10(PDUIPVal[,5]),xlim=c(-.6,.6),ylim=c(0,100), xlab="",ylab="",col=NA,pch=pch0,bg=col0,cex=1.5)
? ? title(xlab="Change in PDUI (Tumor-Normal)",ylab="p-value (BH, -log10)",line= 2.3,cex.lab= 1.3)
? ? abline(v=0.05*c(-1,1),lty=2,lwd=2)
? ? ######################### #cutoff at 0.1改這里#abline(v=0.1*c(-1,1),lty=2,lwd=2)
? ? axis(3,at=c(.55),line=-0.5,label=c("Lengthened"),hadj=1,tick=FALSE)
? ? axis(3,at=c(-.55),line=-0.5, label=c("Shortened"),hadj=0,tick=FALSE)
dev.off()
看下圖：
腫瘤和正常樣本中每個基因的 PDUI 評分圖。虛線代表 ΔPDU=?0.05 截止值粟瞬。藍點代表3'-UTR 延長的基因同仆，紅點代表 3'-UTR 縮短的基因。
?火山圖表示3'-UTR-縮短（紅色）和-延長（藍色）基因（FDR < 0.05）裙品，其 |ΔPDUI|?> 0.05俗批。ΔPDUI （ ΔPDUI =Mean PDUI?T?- Mean PDUI?N）。 |ΔPDUI|?> 0.05 ?作為腫瘤相關(guān)3‘-UTR 縮短或延長的量度（可以更嚴格一點的話市怎，|ΔPDUI|? > 0.1岁忘，但得到的差異基因就更少更少了）。
