對(duì)于集合的可視化，第一時(shí)間想到的都是韋恩圖（venn diagram），一般集合不超過(guò)5個(gè)的時(shí)候，可視化效果還是不錯(cuò)的
但是一旦數(shù)據(jù)集增加，比如說(shuō)五個(gè)的時(shí)候，你就很難從圖中解讀出想要的信息了。
即便你把它畫(huà)的很美觀，如下圖那樣，還是還是很難直觀找到自己需要的信息。可視化的目的不是炫技，而是快速理解數(shù)據(jù)。
還好R語(yǔ)言里新增了一個(gè)集合可視化神包--UpSetR。它可視化的結(jié)果的基礎(chǔ)版本長(zhǎng)下面這個(gè)樣子：
上述是分析了不同電影的所屬類(lèi)型得到的結(jié)果。在我不告訴你任何圖示信息的情況下，請(qǐng)思考下那種電影類(lèi)型拍的最多，然后哪兩種電影電影類(lèi)型拍的最少。
基本上我不用過(guò)多和你解釋圖示，你也能很快的找到答案。圖中黑色表示該位置有數(shù)據(jù)，灰色的點(diǎn)表示沒(méi)有。不同點(diǎn)連線表示存在交集。具體數(shù)據(jù)可以看上面的條形圖。不同類(lèi)型的數(shù)據(jù)的總量看左邊的條形圖。
如何畫(huà)圖
UpSetR是一個(gè)R包，這意味著你可以簡(jiǎn)單通過(guò)一行命令就能安裝
install.packages(UpSetR)
UpsetR接受三種類(lèi)型的數(shù)據(jù)輸入：
- 表格形式，在R語(yǔ)言里就是數(shù)據(jù)框了。行表示元素，列表示數(shù)據(jù)集分配和額外信息。
- 元素名的集合(沒(méi)見(jiàn)過(guò)，不知道。。)
fromList - venneuler包引入的用于描述集合交集的向量
fromExpression。
光看文字肯定是不懂的，所以直接實(shí)戰(zhàn)把
輸入方式一： table
我們用UpSetR提供的測(cè)試數(shù)據(jù)作為演示
require(ggplot2); require(plyr); require(gridExtra); require(grid);
movies <- read.csv(system.file("extdata","movies.csv",package = "UpSetR"), header = TRUE, sep=";")
看下數(shù)據(jù)長(zhǎng)什么樣子
View(movies)
Name是不同的電影，然后不同發(fā)布時(shí)間，后面接著電影跟隨的類(lèi)型。
繪圖用的upset函數(shù)：
upset(movies, nsets = 7, nintersects = 30, mb.ratio = c(0.5, 0.5),
order.by = c("freq", "degree"), decreasing = c(TRUE,FALSE))
稍微解釋一下參數(shù)
nsets: 最多展示多少個(gè)集合數(shù)據(jù)。畢竟原來(lái)有20多種電影類(lèi)型，放不完的
nintersects: 展示多少交集。
mb.ratio： 點(diǎn)點(diǎn)圖和條形圖的比例。
order.by： 交集如何排序。這里先根據(jù)freq，然后根據(jù)degree
decreasing： 變量如何排序。這里表示freq降序，degree升序
更有意思的是，我們還能在圖中描述出1970-1980年恐怖片和劇情片的情況
# 用于query的函數(shù)
between <- function(row, min, max){
newData <- (row["ReleaseDate"] < max) & (row["ReleaseDate"] > min)
}
upset(movies, sets=c("Drama","Comedy","Action","Thriller","Western","Documentary"),
queries = list(list(query = intersects, params = list("Drama", "Thriller")),
list(query = between, params=list(1970,1980), color="red", active=TRUE)))
這里必須介紹一個(gè)神奇的參數(shù)queries：
queries接受query所組成的list。然后不同query也是一個(gè)list，這個(gè)list由查詢(xún)函數(shù)，和參數(shù)組成，參數(shù)也是一個(gè)list。查詢(xún)函數(shù)可以用系統(tǒng)自帶的，也可以自己寫(xiě)一個(gè)。比如說(shuō)這里的between
此外還有一個(gè)參數(shù)叫做attribute.plots能夠添加在upset的結(jié)果圖中加入屬性圖。
upset(movies,attribute.plots=list(gridrows=60,plots=list(list(plot=scatter_plot, x="ReleaseDate", y="AvgRating"),
list(plot=scatter_plot, x="ReleaseDate", y="Watches"),list(plot=scatter_plot, x="Watches", y="AvgRating"),
list(plot=histogram, x="ReleaseDate")), ncols = 2))
這個(gè)attribute.plots接受各個(gè)plot函數(shù)組成的作圖函數(shù)，可以用自帶的，也可以自己寫(xiě)，只要保證里面的參數(shù)設(shè)置正確了。
其他參數(shù)就不繼續(xù)演示了，因?yàn)槲覒小?/p>
輸入方式二：集合交集向量
集合交集向量長(zhǎng)下面這個(gè)樣子
input <- c(
"MAQ"=144600,
"FaSD"=16532,
"Bcftools"=283,
"GATK"=15160,
"MAQ&FaSD"=16323,
"MAQ&Bcftools"=636,
"Bcftools&GATK"=65435,
"FaSD&GATK"=33874,
"MAQ&FaSD&Bcftools"=114,
"MAQ&FaSD&GATK"=41858,
"MAQ&Bcftools&GATK"=4,
"FaSD&Bcftools&GATK"=6603,
"MAQ&FaSD&Bcftools&GATK"=8357
)
輸入格式一目了然，然后數(shù)據(jù)可以用fromExpression進(jìn)行轉(zhuǎn)換
data <- fromExperssion(input)
轉(zhuǎn)換后的數(shù)據(jù)就可以拿去用upset作圖了
upset(data)
福利：Y叔的upsetplot()
我們可以對(duì)ChIP-Seq分析得到的peak進(jìn)行注釋
require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
peakfile <- system.file("extdata", "sample_peaks.txt", package="ChIPseeker")
peakAnno <- annotatePeak(peakfile, tssRegion=c(-3000, 3000), TxDb=txdb)
peakAnno
然后就可以用upsetplot畫(huà)畫(huà)了，太簡(jiǎn)單了。
upsetplot(peakAnno, vennpie=TRUE)
下一期寫(xiě)一篇Y叔的upsetplot是如何寫(xiě)的。
