越来越懒了，不想写记录呀！算了，强迫自己一把！
水平有限，若有错误，恳请指出！
信使RNA ( messenger RNA, mRNA)加工过程中的精细调控对基因的表达具有重要影响,是产生基因功能多样性的重要机制。在前体mRNA成熟过程中,环境的细微变化能够导致在mRNA的不同剪接位点上进行选择性的剪接和多聚腺苷酸化,此过程被称为可选择性多聚腺苷酸化(alerna"tive polyadenylation , APA)。可变多聚腺苷酸化(alternative polyadenylation, APA)是真核细胞mRNA成熟过程中针对前体mRNA 3′端的一种加工修饰方式,是重要的转录后调控机制。
近年来,APA作为调控基因表达的重要方式越来越受到重视。APA现象广泛存在于真核生物中,绝大多数真核基因具有多个多聚腺苷酸[poly( A)]位点。例如,人类有70% ~75%的基因存在APA事件,果蝇1、蠕虫和斑马鱼[3也有约半数的基因存在这类现象。APA能产生长短不同的mRNA转录本,不同的mRNA转录本具有不同的稳定性、翻译能力和表达水平，从而产生多样化的生物学功能。
APA的形成过程APA形成过程包括poly (A)位点的选择和多聚腺苷酸化,这两个过程十分重要并受到多种因素的调控。选择不同的多聚腺苷酸化信号是形成APA的关键。
经典的poly (A)模体是AAUAAA,其他非经典的poly (A)模体变体包括. AUUAAA, AAGAAA和UAUAAA等。在人类基因组中,AAUAAA约占所有多聚腺苷酸化模体的50%,AUUAAA约占16%、A[A/U]UAAA中约20%存在单核苷酸变异,只有约10%的多聚腺苷酸化模体没有类似于AAUAAA的识别区域5mRNA的多聚腺苷酸化包含2个必要步骤:1前体mRNA在poly(A)位点进行剪接;2在3'端加poly (A)尾。
多聚腺苷酸化的主要功能是保证转录本的高效输出、翻译以及稳定。多聚腺苷酸化的顺利形成是前体mRNA序列中的顺式作用元件、反式作用因子以及细胞核中多种酶和蛋白因子相互作用的结果。
选择性多聚腺苷酸化（APA）是一种重要的前体RNA加工机制，广泛存在于所有真核生物中。通过在RNA 3′UTR不同位置上添加polyA尾巴，可以选择性地调节3′UTR的长短。由于3′UTR含有多种顺式调控元件，例如：miRNA或RNA结合蛋白（RBP）结合位点，因此，APA可以通过调控 3′非翻译区(3′UTR) 的长度，影响目标mRNA的稳定性和翻译效率以及翻译后蛋白质的细胞定位，进而精细调节基因表达，对一系列细胞过程（如增殖、分化和肿瘤发生）产生根本性的影响。
内含子多聚腺苷酸化(intronic polyadenylation, IPA)通过形成丢失重要结构域的截短型蛋白实现对靶基因的调控,参与形成肿瘤新抗原。APA具有肿瘤特异性,有可能用于肿瘤分子分型和靶向治疗。
可变多聚腺苷酸化Alternative Polyadenylation (APA)，如下图所示（图片来自参考），在不同的APA信号位点切割，然后添加polyA。这种调控机制属于转录后调控，可能会影响蛋白的序列（发生在编码区），也可能影响蛋白的稳定性（比如非编码区内的miRNA的调控区域）。其实也是可变剪接的一种情况。
一般APA可以分为四种类型：
1 、3’UTR APA：发生在末端外显子内，产生具有不同长度3’UTR的转录本，不影响蛋白编码功能，是最常见的APA形式；
2、 可变末端外显子APA：这种APA产生了末端外显子不同的转录本，影响蛋白编码功能；
3 、内含子APA：发生于在内含子区，延长了某个内部外显子并使之成为末端外显子；
4、 内部外显子APA：在编码区域内部发生剪切和多聚腺苷酸化。
常用的软件是Dapars，这个软件现在也有了升级的版本Dapars2。参考：
https://github.com/ZhengXia/dapars
https://github.com/3UTR/DaPars2
分析流程很相似，Dapars2多了 normalize library sizes 。
安装软件
官网：https://github.com/3UTR/DaPars2
source activate dapars2 # 需要python2.7的环境！
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 脚本 构建注释
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 统计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
### 如果要单独分析每个染色体以节省时间和内存请运行以下命令，可以每条染色体写个脚本，同时提交，同时运行，节省时间和内存 ！
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' 表示转录本的 3'UTR 区域，其他列代表每个样本的 PDUI 值。
再把所有chr 染色体合并到一个文件 all.chr.DaPars2.txt 里。
head all.chr.DaPars2.txt
我有个问题呀，如果我想将 实验组 与 对照组 进行差异分析呢？ 需要怎么实现？我在DaPars2 更新后的版本里，我没有看到这种用法呢？在上一步的configure_file里， 我是将 N 和 T 的wig放在一组运行的。 但怎么做肿瘤组和正常组的差异APA分析吗？ 再研究研究~~
后续来了： 我去给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 后的结果：
在R中，我们取出第一列gene，以及后5列的结果：
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 作为肿瘤相关3‘-UTR 缩短或延长的量度（可以更严格一点的话，|ΔPDUI| > 0.1，但得到的差异基因就更少更少了）。
