CRISPR/Cas9 screen 数据绘图R包
病毒能够侵染细菌，并且巧妙的利用细菌的资源为它所用，可以说是相当的霸道了！简短来说就是病毒就像窃贼一样，能够将自己的基因鸟么悄的整合到细菌基因组上，然后在细菌内大吃二喝疯狂占有资源-利用细菌的元件为自己的基因复制服务。然而细菌也不是逆来顺受的，在其不断的进化过程中建立起一套防盗系统，能够不动声色的把病毒基因从自己的基因组中清除。
细菌的防盗系统拥有多种切除外来病毒基因的功能，科学家们掌握了对一种蛋白Cas9的操作技术，并先后对多种目标细胞DNA进行切除。这种技术被称为CRISPR/Cas9基因编辑系统。今天小编给大家分享的是MAGeCKFlute-哈佛大学刘小乐团队开发的CRISPR下游分析的工具。
简要介绍
MAGeCK (Wei Li and Liu. 2014)和MAGeCK- vispr
刘小乐团队在开发了MAGeCK 和MAGeCK- vispr 用于分析CRISPR/Cas9 screen 数据，为了对目标基因进行下游功能分析，因此开发了MAGeCKFlute。该包的下游分析包括鉴定必要的、非必要的和目标相关的基因，并对这些基因进行生物学功能类别分析、通路富集分析和蛋白复合物富集分析。该包还以多种方式可视化基因，有利于用户探索筛选数据。总的来说，MAGeCKFlute能够准确识别必要的、非必要的和目标基因，以及它们相关的生物学功能。
分析模式
该R包可以进行两种模式的分析，一种是“quick”模式，一种是“step by step”模式。
Quick模式
01
加载R包
library(MAGeCKFlute)
library(ggplot2)02
加载数据
#MAGeCK分析结果数据-gene
file1 = file.path(system.file("extdata", package = "MAGeCKFlute"),
"testdata/rra.gene_summary.txt")
# Read and visualize the file format
gdata = read.delim(file1, check.names = FALSE)
head(gdata)
# id num neg|score neg|p-value neg|fdr neg|rank neg|goodsgrna neg|lfc pos|score pos|p-value pos|fdr pos|rank pos|goodsgrna pos|lfc
#1 CREBBP 4 1.00000 1.00000 1 8269 0 0.96608 3.17e-12 2.9700e-07 0.001650 1 4 0.966082
#2 EP300 4 1.00000 1.00000 1 8270 0 1.02780 3.32e-10 2.9700e-07 0.001650 2 4 1.027803
#3 CHD 4 0.99999 0.99999 1 8268 0 0.59265 1.22e-05 4.4300e-05 0.092203 3 4 0.592654
#4 C16orf72 4 0.99998 0.99998 1 8267 0 0.82307 2.06e-05 8.0000e-05 0.147965 4 4 0.823075
#5 CACNB2 5 0.99863 0.99864 1 8243 0 0.39268 4.45e-05 2.1084e-04 0.319082 5 5 0.392686
#6 FGF12 6 0.54609 0.71111 1 4862 1 0.40822 5.46e-05 2.7627e-04 0.336280 6 5 0.40822
#可以用ReadRRA获取
gdata = ReadRRA(file1)
gdata =ReadRRA(file1, score = c("lfc", "rra")[1])#score 默认选“lfc”，为foldchange值
head(gdata)
# id Score FDR
# 1 CREBBP 0.96608 0.001650
# 2 EP300 1.02780 0.001650
# 3 CHD 0.59265 0.092203
# 4 C16orf72 0.82307 0.147965
# 5 CACNB2 0.39268 0.319082
# 6 FGF12 0.40822 0.336280
#MAGeCK分析结果数据-sgRNA(可选)
file2 = file.path(system.file("extdata", package = "MAGeCKFlute"),
"testdata/rra.sgrna_summary.txt")
sdata = read.delim(file2)
head(sdata)
# sgrna Gene control_count treatment_count control_mean treat_mean LFC control_var adj_var score p.low p.high p.twosided FDR high_in_treatment
#1 s_36798 NF2 8917/21204 5020.7/5127.9 15061.0 5074.30 -1.5693 75491000 78871 35.559 2.9804e-277 1 5.9609e-277 1.2732e-272 False
#2 s_45763 RAB6A 3375.8/3667.7 372.88/357.79 3521.8 365.33 -3.2655 42617 15519 25.338 6.1638e-142 1 1.2328e-141 1.9748e-137 False
#3 s_50164 SF1 3657.8/3352.6 453.62/628.28 3505.2 540.95 -2.6937 46575 15438 23.857 4.2365e-126 1 8.4731e-126 9.0487e-122 False
#4 s_20780 FDXR 3444.8/3191.1 552.51/555.38 3317.9 553.95 -2.5803 32179 14520 22.938 9.7946e-117 1 1.9589e-116 1.5690e-112 False
#5 s_36796 NF2 5492.9/11396 3832.2/3794.6 8444.6 3813.40 -1.1467 17425000 41247 22.803 2.1321e-115 1 4.2641e-115 3.0359e-111 False
#6 s_36799 NF2 1805.2/5541.7 695.86/882.47 3673.5 789.16 -2.2173 6980700 16267 22.615 1.5592e-113 1 3.1183e-113 1.9981e-109 False
sdata = ReadsgRRA(file2，score = c("lfc", "rra")[1])#score 默认选“lfc”，为foldchange值
head(sdata)
# sgrna Gene LFC FDR
# 1 s_36798 NF2 -1.5693 1.2732e-272
# 2 s_45763 RAB6A -3.2655 1.9748e-137
# 3 s_50164 SF1 -2.6937 9.0487e-122
# 4 s_20780 FDXR -2.5803 1.5690e-112
# 5 s_36796 NF2 -1.1467 3.0359e-111
# 6 s_36799 NF2 -2.2173 1.9981e-10903
FluteRRA 分析
#利用FluteRRA对gene和sgRNA数据进行pipeline分析，可接受上边提供的两种文件
FluteRRA(file1, file2, proj="Test", organism="hsa", scale_cutoff = 1, outdir = "./")
# Or
FluteRRA(gdata, sdata, proj="Test", organism="hsa", scale_cutoff = 1, outdir = "./")
#利用FluteRRA仅对gene数据j进行pipeline分析，可接受上边提供的两种文件
FluteRRA(file1, proj="Test", organism="hsa", outdir = "./")
# Or
FluteRRA(gdata, proj="Test", organism="hsa", outdir = "./")
FluteRRA(gdata, proj="Test", organism="hsa", incorporateDepmap = TRUE,
outdir = "./")
Omit common essential genes in the analysis
FluteRRA(gdata, proj="Test", organism="hsa", incorporateDepmap = TRUE,
omitEssential = TRUE, outdir = "./")
#关键参数
#keytype 基因的ID类型，可选"Entrez" or "Symbol"，默认是"Symbol"
#organism="hsa" 可选“hsa”或“mmu”，默认“hsa”
#incorporateDepmap 是否整合Depmap数据，默认“TRUE”
#结果保存在R分析目录下的./MAGeCKFlute_Test/中04
结果展示
./MAGeCKFlute_Test/目录下含有一个pdf文件和RRA文件夹。
1、FluteRRA_Test.pdf 是汇总的结果文件pdf。
2、RRA文件夹下是各个分析的pdf和png以及txt结果文件。
Step by step 模式
01
MAGeCK/MAGeCK-VISPR raw count 结果质控绘图展示
file4 = file.path(system.file("extdata", package = "MAGeCKFlute"),
"testdata/countsummary.txt")
countsummary = read.delim(file4, check.names = FALSE)
head(countsummary)
#File Label Reads Mapped Percentage TotalsgRNAs Zerocounts GiniIndex NegSelQC NegSelQCPval NegSelQCPvalPermutation NegSelQCPvalPermutationFDR NegSelQCGene
#1 ../data/GSC_0131_Day23_Rep1.fastq.gz day23_r1 62818064 39992777 0.6366 64076 57 0.08510 0 1 1 1 0
#2 ../data/GSC_0131_Day0_Rep2.fastq.gz day0_r2 47289074 31709075 0.6705 64076 17 0.07496 0 1 1 1 0
#3 ../data/GSC_0131_Day0_Rep1.fastq.gz day0_r1 51190401 34729858 0.6784 64076 14 0.07335 0 1 1 1 0
#4 ../data/GSC_0131_Day23_Rep2.fastq.gz day23_r2 58686580 37836392 0.6447 64076 51 0.08587 0 1 1 1 0
#Gini index展示
BarView(countsummary, x = "Label", y = "GiniIndex",
ylab = "Gini index", main = "Evenness of sgRNA reads")
# Missed sgRNAs 绘图
countsummary$Missed = log10(countsummary$Zerocounts)
BarView(countsummary, x = "Label", y = "Missed", fill = "#394E80",
ylab = "Log10 missed gRNAs", main = "Missed sgRNAs")
# Read mapping ratio绘图
MapRatesView(countsummary)
02
质控以外的其他分析
这个部分主要包括Depmap相关性分析、过滤一些干扰基因、以及正负选择的可视化。
只是“Quick”模式的拆分，能够对具体的参数进行调整，小编展示一些结果图，具体分析代码可以自己去学习。
正负调控火山图
RRA正负调控排序散点图
功能富集分析图
总结
x阿小编
整体的分析绘图是很方便快速的，两种模式方便选择使用。有兴趣的小伙伴快去测试使用吧！
Binbin Wang, Mei Wang, Wubing Zhang. “Integrative analysis of pooled CRISPR genetic screens using MAGeCKFlute.” Nature Protocols (2019), doi: 10.1038/s41596-018-0113-7.
腾讯云开发者
