bam文件的相关性的比较

bam文件的相关性的比较 by 东林的扯淡小屋

#样本之间的相关性

#根据igv可视化的结果可以初步判断样本之间的相关性很高，所以，我们现在想定量的去计算样本之间的相关性。如何全基因组范围内检查测序数据之间的相关性呢？

#我偏好的方法是按照2000bp划分基因组，然后统计每一个划分模块里面的reads数目，最后做相关性分析，从而定量的检查样本之间的相关性。

cd /data/yudonglin/reference/human/hg19

## 利用samtools统计每一个染色体碱基数目

samtools faidx hg19.fa

## 结果生成一个.fai的文件，其中第一列染色体，第二列是每个染色体碱基的数目

less hg19.fa.fai

awk '{printf $1 "\t" $2 "\n"}' hg19.fa.fai > hg19.genome

## 查看genome文件

head hg19.genome

## 以2000bp为bin划分基因组

bedtools makewindows -g hg19.genome -w 2000 > hg19.windows

## 查看hg19.windows文件，可以看到该文件每行的第二列和第三列之间相差2000，第一列表示染色体

head hg19.windows

## 使用deeptools里面的multiBamSummary功能，统计之前划分的每一个区域在所有样本里面的reads数目

## --numberOfProcessors 4是调用4个线程的意思

multiBamSummary BED-file --BED /data/yudonglin/reference/human/hg19/hg19.windows --bamfiles S1-1_trimmed_sorted.bam S1-2_trimmed_sorted.bam S1-3_trimmed_sorted.bam early1_trimmed_sorted.bam early2_trimmed_sorted.bam early3_trimmed_sorted.bam --numberOfProcessors 40 -out ./readCounts.npz --outRawCounts ./readCounts.tab

## 查看readCounts.tab文件，我们就可以看到三个SRR样本在每一个区间里面的reads数目都被统计出来了，所以接下来我们就可以对这个表格进行相关性分析和PCA分析

less ./readCounts.tab

## 首先下载相关性可视化的包#install.packages("corrplot")## 然后载入这个包library("corrplot")
## 将我们刚才计算得到的数目读进R内（见图28），一共有59837行，6列。第一列是染色体，第二列到第三列是划分的区域，第四列到第五列是之前统计的bam文件名atac_seq <- read.csv("/data1/yudonglin/Rawdata/bowtie2/wang/readCounts.tab",sep = "\t",stringsAsFactors = F)head(atac_seq)## 把前三列删除（这三列不需要，计算相关性只需后面的reads数目就好）atac_seq <- atac_seq[,-1:-3]
## 将剩余的三列的名字重新命名，之前的名字太长了colnames(atac_seq) <- c("S1-1", "S1-2" , "S1-3" , "early1" ,"early2",  "early3")
## 计算样本之间相关性，可以看到样本相关性很强，都在0.99以上（见图29）cor(atac_seq)
## 接着我们使用corrplot包对数据进行可视化，然后将图片以PDF格式输出到本地，然后用AI进行简单的修改就可以放在文章里面了（图30）corrplot(cor(atac_seq))
a<-cor(atac_seq)pheatmap::pheatmap(cor(atac_seq),,scale = "row",                   cluster_rows = T,show_rownames = F,                   cluster_cols = F)

另一种方法：

#样本相似性

conda activate cutandtag

multiBamSummary bins --bamfiles hT-GT_sorted.bam hT-PBS_sorted.bam --binSize 10000 --verbose --extendReads --labels hT-GT hT-PBS -out readCounts.npz --outRawCounts readCounts.tab -p 50

plotCorrelation -in readCounts.npz --removeOutliers --corMethod spearman --skipZeros --whatToPlot heatmap --colorMap jet --plotNumbers -o Corr_readCounts.pdf --outFileCorMatrix Corr_readCounts.tab --plotFileFormat pdf

plotPCA -in readCounts.npz \

-o PCA_readCounts.pdf \

-T "PCA of read counts"

ixxmu / mp_duty

bam文件的相关性的比较 #3577

bam文件的相关性的比较 by 东林的扯淡小屋