学术咖

2021年3月22日2021年4月2日

MetabolAnalyze实现Pareto scaling

BiocManager::install("MetabolAnalyze")
library("MetabolAnalyze")
data <- read.delim('mouse_gc_TMM_cross.xls', row.names = 1, sep = '\t', check.names = FALSE)
col_names <- colnames(data)
row_names <- rownames(data)
data <- as.matrix(data)
b = scaling(data,type = "pareto")
boxplot(b)
write.table(b, file="Pareto_scaled.xls", sep="\t", quote=F, row.names=T, col.names=T)

2021年3月18日2021年3月23日

什么是批次效应（batch effect）？
不同平台的数据，同一平台的不同时期的数据，同一个样品不同试剂的数据，以及同一个样品不同时间的数据等等都会产生一种batch effect 。这种影响如果广泛存在应该被足够重视，否则会导致整个实验和最终的结论失败。比对实验组和对照组，不同的处理是患病和不患病（测序时，先测得疾病，然后测得正常），然后你通过分析，得到很多差异表达的基因。现在问题来了，这个差异表达的结果是和你要研究的因素有关，还是时间有关，这个问题里时间就会成为干扰实验结果的因素，这个效应就是batch effect。

library(limma)
data <- read.table("mouse_gc.xls",sep="\t",header = TRUE)
batch1 <- c(rep('fish',9),rep('mouse',9))
batch1 <- as.factor(batch1)
design <- model.matrix(~0 + batch1)
#data=normalizeBetweenArrays(data)#如果组内的中位数不在同一条水平线上，现在该参数校正，再使用下面的批次校正
new_data <- removeBatchEffect(data[,2:19], batch = batch1)
boxplot(data[,2:19])
boxplot(new_data)
write.table(new_data, file="removeBatch.xls", sep="\t", quote=F, row.names=T, col.names=T)

limma结果出现负值请参考：https://cloud.tencent.com/developer/article/1680305
参考链接：https://www.omicsclass.com/article/1113

2021年3月15日2021年3月20日

preprocessCore分位数标准化

library("preprocessCore")
library("ggplot2")
data <- read.table("mouse_gc.xls",sep="\t",header = TRUE)
rownames(data) <- data$GeneID
#rownames(data) <- data[,1]
#data_mat <- data.matrix(data[,-1]) 
quantile_normalisation <- function(df){
  df_rank <- apply(df,2,rank,ties.method="min")
  df_sorted <- data.frame(apply(df, 2, sort))
  df_mean <- apply(df_sorted, 1, mean)

  index_to_mean <- function(my_index, my_mean){
    return(my_mean[my_index])
  }

  df_final <- apply(df_rank, 2, index_to_mean, my_mean=df_mean)
  rownames(df_final) <- rownames(df)
  return(df_final)
}
new_data <- quantile_normalisation(data[,2:19])
boxplot(data)
boxplot(new_data)
write.table(new_data, file="111.xls", sep="\t", quote=F, row.names=T, col.names=T)
##############################################################
col_names <- colnames(data)
row_names <- rownames(data)
data <- as.matrix(data)
b=normalize.quantiles(data)
boxplot(b)
##使用normalize.quantiles函数和上面结果相同

2021年3月9日2021年4月13日

跨物种比较转录组

近几年几篇高分文献对跨物种比较转录组方法的应用及研究领域，感兴趣的可自行查看原文献：

Cell_2019_跨物种单细胞测序揭示小胶质细胞基因表达谱进化特征

Science_2020_人类猪小鼠大脑中的蛋白编码基因图谱

Science_2018_跨人类和猕猴大脑发育的时空转录组差异

Cell_2021_African lungfish genome sheds lighton the vertebrate water-to-landtransition

Cell_2021_人类大脑进化的发育调控机制

Cell_2021_原始辐鳍鱼类基因组的解析发现登陆的遗传基础已经在硬骨鱼类祖先出现

看完上述文献可知和种内比较转录组相比，跨物种（种间）比较转录组要复杂的多，主要体现在背景基因总数、测序深度、管家基因表达谱和可能的其他因素差异等；分析流程主要分为两步：1、直系同源基因的查找。此处大部分文献使用 one-to-one orthologs （及双向最优比对），如果物种跨度较大可选较为宽松的阈值，期望值1e-10，覆盖度50%；如果物种亲缘关系较为接近就选较严的阈值。也有文献先获得直系同源基因列表，一个gene list及为一个基因（matagene），该基因的表达量为gene list所有基因表达量的和（基因拷贝数差异）。2、合并后基因表达矩阵的标准化。标准化方法各异，方法之一：首先在种内计算出FPKM、RPKM、TPM等，种内TMM标准化，然后根据背景基因集合并数据得到新的表达矩阵再使用TMM标准化。也有直接quantile normalization using the R package preprocessCore（具体R实现可参考下篇文章），标准化方法的选择对结果的影响很大，可自行尝试适合自己的标准化方法。另外，不同数据来源需要去除批次效应，不然进行PCA和样本相关性分析时会将不同数据来源的样本聚到一块；数据的适当缩放也是必要的，一些聚类方法会将极值带来的影响放大化，使结果不准确。下图是Science_2020_人类猪小鼠大脑中的蛋白编码基因图谱分析流程。

2021年2月21日2021年5月16日

数据相关性分析

简介：评价两组数据之间的相关性，有皮尔森（pearson）相关系数，斯皮尔曼（spearman）相关系数和肯德尔（kendall）相关系数。在这三大相关系数中，spearman和kendall属于等级相关系数亦称为“秩相关系数”，是反映等级相关程度的统计分析指标。相关系数的绝对值越大，相关性越强，相关系数越接近于1或-1，相关度越强，相关系数越接近于0，相关度越弱。

Pearson（皮尔逊）相关系数：皮尔逊相关也称为积差相关（或积矩相关）是英国统计学家皮尔逊于20世纪提出的一种计算直线相关的方法。适用于：

(1)、两个变量之间是线性关系，都是连续数据。

(2)、两个变量的总体是正态分布，或接近正态的单峰分布。

(3)、两个变量的观测值是成对的，每对观测值之间相互独立。

Spearman Rank（斯皮尔曼等级）相关系数：在统计学中，斯皮尔曼等级相关系数以Charles Spearman命名，并经常用希腊字母ρ（rho）表示其值。又称秩相关系数，是利用两变量的秩次大小作线性相关分析，对原始变量的分布不作要求，属于非参数统计方法，适用范围要广些。斯皮尔曼等级相关是根据等级资料研究两个变量间相关关系的方法。它是依据两列成对等级的各对等级数之差来进行计算的，所以又称为“等级差数法”。斯皮尔曼等级相关对数据条件的要求没有积差相关系数严格，只要两个变量的观测值是成对的等级评定资料，或者是由连续变量观测资料转化得到的等级资料，不论两个变量的总体分布形态、样本容量的大小如何，都可以用斯皮尔曼等级相关来进行研究。对于服从Pearson 相关系数的数据亦可计算Spearman 相关系数，但统计效能要低一些。Pearson 相关系数的计算公式可以完全套用 Spearman 相关系数计算公式，但公式中的x 和y 用相应的秩次代替即可。

Kendall Rank（肯德尔等级）相关系数：在统计学中，肯德尔相关系数是以Maurice Kendall命名的，并经常用希腊字母τ（tau）表示其值。用于反映分类变量相关性的指标，适用于两个分类变量均为有序分类的情况。对相关的有序变量进行非参数相关检验；取值范围在-1-1 之间，此检验适合于正方形表格；肯德尔(Kendall)W 系数又称和谐系数，是表示多列等级变量相关程度的一种方法。适用这种方法的数据资料一般是采用等级评定的方法收集的，即让K 个评委（被试）评定N 件事物，或1 个评委（被试）先后K 次评定N 件事物。等级评定法每个评价者对N 件事物排出一个等级顺序，最小的等级序数为1 ，最大的为N ，若并列等级时，则平分共同应该占据的等级，如，平时所说的两个并列第一名，他们应该占据1 ，2 名，所以它们的等级应是1.5, 又如一个第一名，两个并列第二名，三个并列第三名，则它们对应的等级应该是1,2.5,2.5,5,5,5, 这里2.5 是2,3 的平均，5 是4,5,6 的平均。

R包计算相关性系数：

library(ggplot2)
library(ggpubr)
data <- read.table("111.xls",sep="\t",header = TRUE)#列是特征（两个特征），行是物种名
ggplot(data=data, aes(x=genome, y=Tes))+geom_point(color="red")+stat_smooth(method="lm")+stat_cor(data=data, method = "pearson")+theme_classic()
#stat_cor(data=dat, method = "pearson")意为用pearson相关进行相关性分析，可以自行更改方法
#stat_smooth是画拟合曲线的函数
#se=FALSE意思为不画出置信区间
#data=后跟需要画图的数据的文件名
#X=后跟作为X轴的数据的那一列的列名
#Y=后跟作为Y轴的数据的那一列的列名
#geom_point函数是个性化设置散点图点的形状，颜色，大小等，此处只设置了颜色，有需要可自行加入
#theme_classic() x,y轴实线化
library(corrplot)
library(pheatmap)
corr <- cor(data, method = 'pearson')
pheatmap(corr)
corrplot(corr, type = 'upper', tl.col = 'black', order = 'hclust', tl.srt = 45, addCoef.col = 'white')

2021年2月20日2021年2月25日

MUMmer使用及后续绘图

MUMmer被广泛用于大片段序列的比对，如染色体共线性分析。

nucmer [options]  
delta-filter -i 80 -l 1000 -r -q out.delta > out.rq.delta
#比对率大于80%，对比长度大于1000
#-r: 仅保留每个reference在query上的最佳位置,允许多条reference在query上重叠
#-q: 仅保留每个query在reference上的最佳位置,允许多条query在reference上重叠
show-coords out.rq.delta > out.coords
grep -P  "zfCh04\s+" out.coords|awk '{print $12,$13}' |sort |uniq -c #查看reference单个染色体zfCh04和query不同染色体的匹配区域的数量
grep -P  "zfCh04\s+gcCh04" out.coords|awk '{print $12,$1,$2,$13,$4,$5}' > zfCh04_gcCh04.txt
sed -i 's/ /\t/g' zfCh04_gcCh04.txt

#R作图
使用RIdeogram包，可参考RIdeogram：染色体数据可视化的R包

install.packages('RIdeogram')
require(RIdeogram)
cc <- read.table("111.xls",sep="\t",header = TRUE,stringsAsFactors = F)
dd <- read.table("222.xls",sep="\t",header = TRUE,stringsAsFactors = F)
ideogram(karyotype = dd, synteny = cc)
#data(karyotype_dual_comparison, package="RIdeogram")
#data(synteny_dual_comparison, package="RIdeogram")
ideogram(karyotype = karyotype_dual_comparison, synteny = synteny_dual_comparison)
convertSVG("chromosome.svg", device = "png")

karyotype_dual_comparison文件格式
Chr: 染色体号
Start: 起始
End: 终止
fill: 染色体填充色
species：物种名
size: 物种名字体大小
color: 物种名字体颜色

synteny_dual_comparison文件格式
Species_1：物种1染色体号
Start_1，End_1：物种1染色体区域位置
Species_2：物种2染色体号
Start_2，End_2：物种2染色体区域位置

此外还支持三个基因组的共线性

2021年1月6日2021年1月6日

ATAC-seq分析流程（三）

前边介绍了使用DiffBind进行ATAC-seq数据的差异分析，DiffBind进行差异分析时调用了edgeR或者DESeq2包；相对于DiffBind流程我更习惯普通转录组分析流程，因为两者本质上是一样的，区别在于ATAC-seq要根据结果的peak文件自己制备gtf文件，peak对应的gtf文件和gene是不一样的，不同样本间也可能存在差异。

cat *.narrowPeak > merged.peaks
sort -k1,1 -k2,2n merged.peaks > merged.peaks.sorted
bedtools merge -i merged.peaks.sorted > bedtools.merged.sorted.bed
awk -v var=ATAC_seq '{print $1"\t"var"\texon\t"$2"\t"$3"\t.\t+\t.\tgene_id \"peak_"NR"\"; transcript_id \"peak_"NR"\";"}' bedtools.merged.sorted.bed > bedtools.merged.atac.gtf
featureCounts -T 18 -p -t exon -g gene_id -a bedtools.merged.atac.gtf -o ATAC_counts.txt *last.bam
sed 1d ATAC_counts.txt | cut -f 1,7-15 > rawATAC_counts.matrix

一个peak相当于一个gene，得到的原始矩阵可以使用DESeq2包分析，筛选出共有特有或者差异的peaks等进行peaks的注释。

2021年1月5日2021年4月9日

ATAC-seq分析流程（二）

六、Deeptools可视化

#deeptools安装conda install deeptools
gtf2bed < ~/gc_genome/17.Geta/1.geta/out.gtf > genome.bed
for i in `cat ../samples.txt`; do echo "bamCoverage -p 5 --normalizeUsing RPKM -b $i.last.bam -o $i.last.bw"; done > bam2bw.list
ParaFly -c bam2bw.list -CPU 9
computeMatrix reference-point  --referencePoint TSS  -p 15 -b 10000 -a 10000 -R ../genome.bed -S ../*bw --skipZeros  -o matrix1_test_TSS.gz --outFileSortedRegions regions1_test_genes.bed
plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.pdf --plotFileFormat pdf  --dpi 720
plotProfile -m matrix1_test_TSS.gz  -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720
computeMatrix scale-regions -p 150 -b 10000 -a 10000 -R ../genome.bed -S ../*bw --skipZeros  -o matrix1_test_body.gz --outFileSortedRegions regions1_test_body.bed
plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.pdf --plotFileFormat pdf  --dpi 720
plotProfile -m matrix1_test_TSS.gz  -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720

七、重复样品处理（至少两个生物学重复）
IDR是通过比较一对经过排序的regions/peaks 的列表，然后计算反映其重复性的值，同时得到了merge后的一致性的peaks。

for i in `cat ../../samples.txt`; do echo "sort -k8,8nr ../$i.peaks_peaks.narrowPeak > $i.peaks.narrowPeak"; done > peaks_sort.list
ParaFly -c peaks_sort.list -CPU 9

八、peaks注释（非模式物种）
使用GenomicFeatures包的函数制作TxDb对象：

#makeTxDbFromUCSC： 通过UCSC在线制作TxDb
#makeTxDbFromBiomart: 通过ensembl在线制作TxDb
#makeTxDbFromGRanges：通过GRanges对象制作TxDb
#makeTxDbFromGFF：通过解析GFF文件制作TxDb
gc_TxDB <- makeTxDbFromGFF("genome.gff3")
library(clusterProfiler)
library("ChIPseeker")
library(ChIPpeakAnno)
library(GenomicFeatures)
library("ggupset")
library("ggimage")
gc_TxDB <- makeTxDbFromGFF("out.gff3")
promoter <- getPromoters(TxDb=gc_TxDB, upstream=3000, downstream=3000)
files = list(A_summits = "A.peaks_summits.bed", B_summits = "B.peaks_summits.bed", C_summits = "C.peaks_summits.bed", D_summits = "D.peaks_summits.bed")
peakAnno <- annotatePeak(files[[1]], tssRegion=c(-3000, 3000), TxDb=gc_TxDB)
plotAnnoPie(peakAnno)
#tiff("peakAnno.tiff")
#plotAnnoPie(peakAnno)
#dev.off() 此方法保存的图片标注有颜色模块
plotAnnoBar(peakAnno)
vennpie(peakAnno)
peakAnnoList <- lapply(files, annotatePeak, TxDb=gc_TxDB,tssRegion=c(-3000, 3000))
plotAnnoBar(peakAnnoList)
#tagHeatmap来画热图
peak <- readPeakFile(files[[3]])
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
peakHeatmap(files, TxDb=gc_TxDB, 
    upstream=3000, downstream=3000, 
    color=rainbow(length(files)))
#以谱图的形式来展示结合的强度
plotAvgProf2(files[[3]], TxDb=gc_TxDB, 
			 upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency")
tagMatrixList <- lapply(files, getTagMatrix, 
                        windows=promoter)
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))
#注释基因
peakAnnoList <- lapply(files, annotatePeak, TxDb=gc_TxDB,tssRegion=c(-3000, 3000),addFlankGeneInfo=TRUE, flankDistance=5000)
genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)
save(peakAnnoList,file="peakAnnolist.rda")
write.table(as.data.frame(peakAnnoList$mulvscon),file="Nanog.PeakAnno",sep='\t',quote = F)

九、差异分析
DiffBind安装

conda install libv8
conda install -c conda-forge librsvg
##R version 4.0.3 (2020-10-10)
BiocManager::install("systemPipeR")
BiocManager::install("DiffBind")
library(DiffBind)
setwd("D:/Experiment_data/wcs/ATAC/")
tamoxifen <- dba(sampleSheet="sheet_ATAC.CSV")
tamoxifen <- dba.count(tamoxifen, bUseSummarizeOverlaps=TRUE)
plot(tamoxifen)
tamoxifen <- dba.contrast(tamoxifen, categories = DBA_FACTOR, minMembers = 3)
tamoxifen <- dba.analyze(tamoxifen, method = DBA_DESEQ2)
#heatmap
hmap <- colorRampPalette(c("green", "black", "red"))(n = 13)
dba.plotHeatmap(tamoxifen, correlations=FALSE, scale="row", colScheme = hmap)
#MA plots
dba.plotMA(tamoxifen, contrast = 2)
#火山图
dba.plotVolcano(tamoxifen)
tamoxifen.DB <- dba.report(tamoxifen,  method=DBA_DESEQ2, contrast = 1, th=1)
#write.csv(tamoxifen.DB,"diffbind.csv")
out <- as.data.frame(tamoxifen.DB)
write.table(out, file="tamoxifen_DB.csv", sep=",", quote=F, col.names = NA)
#保存不同样本的表达矩阵
out <- as.data.frame(tamoxifen[["binding"]])
write.table(out, file="deseq2_count.xls", sep="\t", quote=F, row.names=F, col.names=F)
#以bed格式保存显著性的差异结果
out <- as.data.frame(tamoxifen.DB)
deseq.bed <- out[ which(out$FDR < 0.05), 
                 c("seqnames", "start", "end", "strand", "Fold")]
write.table(deseq.bed, file="results/contrast_1_deseq2_sig.bed", sep="\t", quote=F, row.names=F, col.names=F)

参考链接：https://github.com/ENCODE-DCC/atac-seq-pipeline

文章

Gene-rank score计算