在线生信绘图-各种生信绘图软件的区别

大多数科研文章都离不开图表,尤其是图,熟悉一些绘图软件,并将图在文章和PPT中展示出来,是科研训练的重要内容。漂亮的文章配图能给自己的工作加不少分,生信宝典推出R的系列教程 ggplot2高效实用指南 (可视化脚本、工具、套路、配色) 讲解通过R语言绘制高颜值图。后来为了更加方便使用,生信宝典团队开发了在线绘图工具 ,支持14中常见图形和部分扩增子分析,深受欢迎,日均访问400次,累计访问数十万次,遍及世界各大洲,功能也在一直增加完善中。

后台不少朋友说有更多绘图的需求,故整理了

Graphpad

的使用,一是传播知识,二是向优秀软件学习,让 在线绘图 越来越好。

Graphpad , Origin和SigmaPlot

是科研领常用的三款统计作图软件,他们各有优缺点,经常需要配合使用。

今天给大家介绍的是Graphpad的基本使用方法

下面是Graphpad的欢迎界面

从

方框1

处可知Graphpad Prism 6为用户提供了

6种

基本统计图表选择;

方框2

处是对所选类型图表的解释说明,点击 Learn more 查看官方帮助文档;

方框3

处可以选择某一类型图表下的小分类;

方框4

处则提供了演示数据集,方便新手入门。

6种统计图表类型的说明和使用的统计方法介绍如下:

在此选择 grouped 下提供的第二个现成数据集进行作图演示。选好后直接点击 creat 即可。

上一步之后我们便进入了graphpad操作的主界面,

方框1

处展示的是我们的数据,每个大组表示三种不同的细胞系(wild-type cells、GPP5 cell line和GPP7 cell line);每行表示不同的处理(血清饥饿和正常培养);每组5个平行重复。

该实验数据存在缺失值,如第一行的A:Y4。

方框2

处是一个

目录树

导入数据后,点击 Graphs 即可出图,如下:

在这个界面中,可以根据自己的需求改变图的类型。初学者要敢于

乱点

,自己试试不同的选项,就能知道每部分是什么功能了。( 尝试出真知 )

上面初步作出的统计图并没有展示诸如组间差异的统计分析结果,那我们要怎么添加 * 以展示组间差异呢?

在欢迎界面有 enter and plot error values already calculated elsewhere 选项,说明可以从外部导入在别处做好的统计分析结果,那是要去其他地方做好计算再来作图吗?

不需要!因为Graphpad上便可直接做统计分析!!!虽然可以在别的地方先做好分析,但是并不建议这样做,直接一个操作更便捷。

首先确定我们的分析目标:

1)评估细胞系之间的差异是否超过预期;

2)评估处理之间的差异是否超过预期的偶然性;

3)评估处理之间的差异是否与每种细胞系一致;

4)计算不同处理的细胞系之间差异的95%置信区间。

因此我们需要进行

两因素的方差分析

,点击 Analyze ,从Grouped分析列表中选择 Two-way ANOVA ,并接受对话框中的所有默认选项。Graphpad人性化之处在于,会自动识别当前分析的数据结构,给用户展示出最优的默认分析选项。

之后的具体分析便需要一定的统计学知识了,初生牛犊别怕试,你会发现点击不同的选项时,图示也会跟着变化,十分方便用户理解。

1)第一部分:实验设计,最开始提到了该组数据是非配对的,选第一;

2)第二部分:多重比较选择,就是你想要横着比还是竖着比或是斜着比;

3)第三部分:选择做统计分析的一些参数。

上面都设置好后,点击 OK ,分析结果会展示到操作主界面,并且保存到左侧 目录树 下的 Results 文件夹中。

1. 图形编辑工具栏介绍

1)第一部分:文档操作模块,跟office类似,如打开、新建、保存、复制粘贴等;

2)第二部分:数据分析,即之前演示的对录入数据进行统计分析;

3)第三部分:图形编辑,这部分是该软件的精华所在。可对图形进行例如图形类型、坐标轴、误差棒、图形大小、颜色等编辑更改;

4)第四部分:文字编辑,可以在统计图上插图线条、添加文字等;

5)第五部分:图形输出,可将做好的图导出、打印或输出到word/PPT中。

2. 根据之前的统计分析结果,利用文字编辑工具给统计图添上组间差异的标志

1)draw:可以插入不同形状的线;

2)write:编辑文字,插入公式等;

3)双击图形中的任一元素,可对文字、线条等进行形状、大小、粗细等编辑

3. 坐标轴调整

点击 Change 中的

第二个图标

,或者在图中

双击坐标轴

,调出坐标调整工具;

坐标轴整体设置和X轴设置,见图解,一般不做太多改动;

重点讲解如何设置Y轴,做出

截断图

,这是统计作图中常见的问题。

之前的数据相差不大,不太适合做演示,故另外生成了一个柱状图。左边柱子值太大,右边柱子值太小,右侧的柱子几乎消失了,使图看起来很不协调。做一个

截断图

可以很好地解决该问题。

1)点击坐标设置中的 Left Y axis (一般作图默认只生成左边的Y轴);

2) Gaps and Direction 处选择 Two segments ,即将Y轴断裂成两部分(根据实际数据,还可以选择将Y轴断裂成3段);

3)在 Segment 处,对下段(Bottom)和上段(Top)的Y轴分别设置参数(Rang);

4)然后在 Regularly spaced ticks 下的 Major ticks interval 处设置Y轴数值标签;

length:50%,上下段的长度各一半;

minimum:下段最小值0,上段最小值150;

maximum:下段最大值10,上段最大值200;

Major ticks:下段每格大小为2,上段每格大小为10。

这样便将Y轴10-150中间部分截掉了,整个图便显得好看多了。

4. 颜色修改

点击 Change 中的

最后一个图标

,或者在图中

双击柱子

,调出调整参数的工具框;

利用graphpad做统计图的基本操作大概就是以上这些内容,掌握一个工具的使用方法并不难,多试试便会了。难的是背后的统计学知识,如何对自己的数据选择合适的统计分析方法?当掌握不同分析方法适用于哪些数据结构后,便可以举一反三,对不同实验设计的数据做合适的分析了。

各种生信绘图软件的区别就是软件的用户群体不同,对于不同的用户,软件会对用户的行为进行改变。

ATAC-seq信息分析流程主要分为以下几个部分:数据质控、序列比对、峰检测、motif分析、峰注释、富集分析,下面将对各部分内容进行展开讲解。

下机数据经过过滤去除接头含量过高或低质量的reads,得到clean reads用于后续分析。常见的trim软件有Trimmomatic、Skewer、fastp等。fastp是一款比较新的软件,使用时可以用--adapter_sequence/--adapter_sequence_r2参数传入接头序列,也可以不填这两个参数,软件会自动识别接头并进行剪切。如:

fastp \

--in1 A1_1. \ # read1原始fq文件

--out1 A1_clean_1. \ # read1过滤后输出的fq文件

--in2 A1_2.  \ # read2原始fq文件

--out2 A1_clean_2. \ # read2过滤后输出的fq文件

--cut_tail  \ #从3’端向5’端滑窗,如果窗口内碱基的平均质量值小于设定阈值,则剪切

--cut_tail_window_size=1 \ #窗口大小

--cut_tail_mean_quality=30 \ #cut_tail参数对应的平均质量阈值

--average_qual=30 \ #如果一条read的碱基平均质量值小于该值即会被舍弃

--length_required=20  \ #经过剪切后的reads长度如果小于该值会被舍弃

fastp软件的详细使用方法可参考:

。fastp软件对于trim结果会生成网页版的报告,可参考官网示例

和

,也可以用FastQC软件对trim前后的数据质量进行评估,FastQC软件会对单端的数据给出结果,如果是PE测序需要分别运行两次来评估read1和read2的数据质量。

如:

fastqc A1_1.

fastqc A1_2.

FastQC会对reads从碱基质量、接头含量、N含量、高重复序列等多个方面对reads质量进行评估,生成详细的网页版报告,可参考官网示例:

经过trim得到的reads可以使用BWA、bowtie2等软件进行比对。首先需要确定参考基因组fa文件,对fa文件建立索引。不同的软件有各自建立索引的命令,BWA软件可以参考如下方式建立索引:

bwa index

建立好索引后即可开始比对,ATAC-seq推荐使用mem算法,输出文件经samtools排序输出bam:

bwa mem  A1_clean_1. A1_clean_2.

| samtools sort -O bam -T A1 >

值得注意的是,在实验过程中质体并不能完全去除,因此会有部分reads比对到质体序列上,需要去除比对到质体上的序列,去除质体序列可以通过samtools提取,具体方法如下:首先将不含质体的染色体名称写到一个chrlist文件中,一条染色体的名称写成一行,然后执行如下命令即可得到去除质体的bam

samtools view -b $chrlist > _MT_

用于后续分析的reads需要时唯一比对且去重复的,bwa比对结果可以通过MAPQ值来提取唯一比对reads,可以用picard、sambamba等软件去除dup,最终得到唯一比对且去重复的bam文件。

比对后得到的bam文件可以转化为bigWig(bw)格式,通过可视化软件进行展示。deeptools软件可以实现bw格式转化和可视化展示。首先需要在linux环境中安装deeptools软件,可以用以下命令实现bam向bw格式的转换:

bamCoverage -b -o

此外,可以使用deeptools软件展示reads在特定区域的分布,如:

computeMatrix reference-point   \ # reference-pioint表示计算一个参照点附近的reads分布,与之相对的是scale-regions,计算一个区域附近的reads分布

--referencePoint TSS   \#以输入的bed文件的起始位置作为参照点

-S   \ #可以是一个或多个bw文件

-R   \ #基因组位置文件

-b 3000   \ #计算边界为参考点上游3000bp

-a 3000   \ #计算边界为参考点下游3000bp,与-b合起来就是绘制参考点上下游3000bp以内的reads分布

-o   \ #输出作图数据名称

#图形绘制

plotHeatmap \

-m  new_ \ #上一步生成的作图数据

-out \ # 输出图片名称

绘图结果展示:

MACS2能够检测DNA片断的富集区域,是ATAC-seq数据call peak的主流软件。峰检出的原理如下:首先将所有的reads都向3'方向延伸插入片段长度,然后将基因组进行滑窗,计算该窗口的dynamic λ,λ的计算公式为:λlocal = λBG(λBG是指背景区域上的reads数目),然后利用泊松分布模型的公式计算该窗口的显著性P值,最后对每一个窗口的显著性P值进行FDR校正。默认校正后的P值(即qvalue)小于或者等于0.05的区域为peak区域。需要现在linux环境中安装macs2软件,然后执行以下命令:

macs2 callpeak \

-t \ #bam文件

-n A1 \ # 输出文件前缀名

--shift -100 \ #extsize的一半乘以-1

--extsize 200 \ #一般是核小体大小

--call-summits #检测峰顶信息

注:以上参数参考文献(Jie Wang,.2018.“ATAC-Seq analysis reveals a widespread decrease of chromatin accessibility in age-related macular degeneration.”Nature Communications)

ATAC分析得到的peak是染色质上的开放区域,这些染色质开放区域常常预示着转录因子的结合,因此对peak区域进行motif分析很有意义。常见的motif分析软件有homer和MEME。以homer软件为例,首先在linux环境中安装homer,然后用以下命令进行motif分析:

\

A1_ \ #用于进行motif分析的bed文件

 \ #参考基因组fa文件

A1  \ #输出文件前缀

-size  given \ #使用给定的bed区域位置进行分析,如果填-size -100,50则是用给定bed中间位置的上游100bp到下游50bp的区域进行分析

homer分析motif的原理及结果参见:

根据motif与已知转录因子的富集情况可以绘制气泡图,从而可以看到样本与已知转录因子的富集显著性。

差异peak代表着比较组合染色质开放性有差异的位点,ChIP-seq和ATAC-seq都可以用DiffBind进行差异分析。DiffBind通过可以通过bam文件和peak的bed文件计算出peak区域标准化的readcount,可以选择edgeR、DESeq2等模型进行差异分析。

在科研分析中我们往往需要将peak区域与基因联系起来,也就是通过对peak进行注释找到peak相关基因。常见的peak注释软件有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker为例,需要在R中安装ChIPseeker包和GenomicFeatures包,然后就可以进行分析了。

library(ChIPseeker)

library(GenomicFeatures)

txdb<- makeTxDbFromGFF(‘’)#生成txdb对象,如果研究物种没有已知的TxDb,可以用GenomicFeatures中的函数生成

peakfile <-readPeakFile(‘A1_’)#导入需要注释的peak文件

peakAnno <- annotatePeak(peakfile,tssRegion=c(-2000, 2000), TxDb=txdb)

# 用peak文件和txdb进行peak注释,这里可以通过tssRegion定义TSS区域的区间

对于peak注释的结果,也可以进行可视化展示,如:

p <- plotAnnoPie(peakAnno)

通过注释得到的peak相关基因可以使用goseq、topGO等R包进行GO富集分析,用kobas进行kegg富集分析,也可以使用DAVID在线工具来完成富集分析。可以通过挑选感兴趣的GO term或pathway进一步筛选候选基因。