精华内容
下载资源
问答
  • 通过降维技术把多个变量化为少数几个主成分(综合变量)的统计分析方法。这些主成分能够反映原始变量的绝大部分信息,它们通常表示为原始变量的某种线性组合。 主成分分析经常用于减少数据集的维数,同时保持数据...

     

    主成分分析(英语:Principal components analysis,PCA)是一种分析、简化数据集的技术。

    通过降维技术把多个变量化为少数几个主成分(综合变量)的统计分析方法。这些主成分能够反映原始变量的绝大部分信息,它们通常表示为原始变量的某种线性组合。 

    主成分分析经常用于减少数据集的维数,同时保持数据集中的对方差贡献最大的特征。这是通过保留低阶主成分,忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。

    主成分分析由卡尔•皮尔逊于1901年发明,用于分析数据及建立数理模型。其方法主要是通过对协方差矩阵进行特征分解,以得出数据的主成分(即特征向量)与它们的权值(即特征值)。

    主成分的目的: 


    (1)变量的降维 
    (2)主成分的解释(在主成分有意义的情况下)

    主成分分析法从冗余特征中提取主要成分,在不太损失模型质量的情况下,提升了模型训练速度。

    如上图所示,我们将样本到红色向量的距离称作是投影误差(Projection Error)。以二维投影到一维为例,PCA 就是要找寻一条直线,使得各个特征的投影误差足够小,这样才能尽可能的保留原特征具有的信息。因为PCA仅保留了特征的主成分,所以PCA是一种有损的压缩方式.

    PCA分析的一般步骤

    1.根据研究问题选取初始分析变量

    2.根据初始变量特性判断由协方差阵求主成分还是由相关矩阵求主成分;

    3.求协方差阵或相关阵的特征值与相应标准特征向量;

    4.判断是否存在明显的多重共线性,若存在,则回到第(1)步;

    5.得到主成分的表达式并确定主成分个数,选取主成分;

    6.结合主成分对研究问题进行分析并深入研究。

    PCA的目标是用一组较少的不相关变量代替大量相关变量,同时尽可能保留初始变量的信息,这些推导所得的变量称为主成分,它们是观测变量的线性组合。

    主成分分析法优缺点

    优点

    ↘可消除评估指标之间的相关影响。因为主成分分析法在对原始数据指标变量进行变换后形成了彼此相互独立的主成分,而且实践证明指标间相关程度越高,主成分分析效果越好。

    ↘可减少指标选择的工作量,对于其他评估方法,由于难以消除评估指标间的相关影响,所以选择指标时要花费不少精力,而主成分分析法由于可以消除这种相关影响,所以在指标选择上相对容易些。

    ↘主成分分析中各主成分是按方差大小依次排列顺序的,在分析问题时,可以舍弃一部分主成分,只取前面方差较大的几个主成分来代表原变量,从而减少了计算工作量。用主成分分析法作综合评估时,由于选择的原则是累计贡献率≥85%,不至于因为节省了工作量却把关键指标漏掉而影响评估结果。

    缺点

    ↘在主成分分析中,我们首先应保证所提取的前几个主成分的累计贡献率达到一个较高的水平(即变量降维后的信息量须保持在一个较高水平上),其次对这些被提取的主成分必须都能够给出符合实际背景和意义的解释(否则主成分将空有信息量而无实际含义)。

    ↘主成分的解释其含义一般多少带有点模糊性,不像原始变量的含义那么清楚、确切,这是变量降维过程中不得不付出的代价。因此,提取的主成分个数m通常应明显小于原始变量个数p(除非p本身较小),否则维数降低的“利”可能抵不过主成分含义不如原始变量清楚的“弊”。

    ↘当主成分的因子负荷的符号有正有负时,综合评价函数意义就不明确。

    主成分分析案例

    某公司经理拟招聘一名员工,要求其具有较高的工作积极性、自主性、热情和责任感。为此,该经理专门设计了一个测试问卷,配有25项相关问题,拟从315位应聘者中寻找出最合适的候选人。

    在这25项相关问题中:

    ↘Qu3-Qu8、Qu12、Qu13测量的是工作积极性

    ↘Qu2、Qu14-Qu19测量工作自主性

    ↘Qu20-Qu25测量的是工作热情

    ↘Qu1、Qu9-Qu11测量工作责任感

    每一个问题都有非常同意“Agree”、同意 “Agree Some”、不确定“Undecided”、不同意 “Disagree Some”和 非常不同意 “Disagree”五个等级。

    该经理想根据这25项问题判断应聘者在这四个方面的能力,现收集了应聘者的问卷信息,经汇总整理后部分数据如下:

    分析者希望将多个变量归纳为某几项信息进行分析,即降低数据结果的维度。针对这种情况,可以进行主成分提取,但需要先满足2项假设:

    ↘假设1:观测变量是连续变量或有序分类变量,如本研究中的测量变量都是有序分类变量。

    ↘假设2:变量之间存在线性相关关系。

    SPSS操作

    SPSS操作

    (1) 在主页面点击Analyze→Dimension Reduction →Factor

    (2) 将变量Qu1-Qu25放入Variables栏

    (3) 点击Deive,点选Statistics栏的Initial solution选项,并点选Correlation Matrix栏的Coefficients、KMO and Bartlett’s test of sphericity、Reproduced和Anti_image选项

    (4) 点击Continue→Extraction,点击Display栏中的Scree plot选项

    (5) 点击Continue→Rotation,点选Method栏的Varimax选项,并点选Display栏的Rotated solution和Loading plot(s)选项

    (6) 点击Continue→Scores,点击Save as variables,激活Method栏后点击Regression选项

    (7) 点击Continue→Options,点击 Sorted by size和Suppress small coefficients选项,在Absolute value below栏内输入“.3”点击Continue→OK

    经上述操作,SPSS输出相关矩阵表如下:

    该表主要用于判断各变量之间的线性相关关系,从而决定变量的取舍,即如果某一个变量与同一分组中其他变量之间的关联性不强,我们就认为该变量与其他变量测量的内容不同,在主成分提取中不应该纳入该变量。一般来说,如果相关系数大于等于0.3,我们就认为变量之间存在较好的线性相关性。

    从本研究的结果来看,在分别对应聘者工作积极性(Q3-Q8,Q12,Q13)、工作自主性 (Q2,Q14-19)、工作热情(Q20-25)和工作责任感(Q1,Q9-11)的测量中,每组变量之间的相关系数均大于0.3,说明各组变量之间具有线性相关关系,提示满足假设2。

    KMO检验对数据结构的总体分析

    KMO检验主要用于主成分提取的数据情况。KMO检验系数分布在0到1之间,如果系数值大于0.6,则认为样本符合数据结构合理的要求。

    部分学者认为,只有当KMO检验系数值大于0.8时,主成分分析的结果才具有较好的实用性,具体系数对应关系如下:

    SPSS输出本研究结果如下:

    本研究的KMO检验系数为0.833,根据系数对应关系表,我们认为本研究数据结构很好(meritorious),具有相关关系,满足假设2。

    KMO检验对各变量的单独分析

    SPSS输出各变量的KMO检验结果如下:

    整理后各题KMO值:

    KMO检验对单个变量的分析结果也在0到1之间分布,如果系数大于0.5,则认为单个变量满足要求;如果系数大于0.8,则认为单个变量结果很好。

    分析结论中,任一变量的KMO检验结果均大于0.7,即各变量结果一般,但满足假设2。

    Bartlett's检验

    Bartlett's检

    Bartlett's检验的零假设是研究数据之间的相关矩阵是一个完美矩阵,即所有对角线上的系数为1,非对角线上的系数均为0。

    在完美矩阵情况下,各变量之间没有相关关系,即不能将多个变量简化为少数的成分,没有进行主成分提取的必要。因此,我们希望拒绝Bartlett's检验的零假设。

    SPSS输出结果如下:

    Bartlett's检验的P值小于0.001,拒绝零假设,即认为研究数据可以进行主成分提取,满足假设2。

    结果解释

    对主成分结果的分析主要从公因子方差(communalities)、提取主成分和强制提取主成分三个方面进行。

    公因子方差结果

    SPSS输出公因子方差结果如下:

    研究中有多少个变量数据结果就会输出多少个成分,本研究中共有25个变量,就会对应产生25个成分。

    “Extraction”栏提示当只保留选中的成分时,变量变异被解释的程度。

    提取主成分

    研究中有多少个变量,主成分提取就会产生多少个主成分。我们通过选取主成分对数据进行降维,但同时也要注意尽可能多地包含对数据变异的解释。

    一般来说,结果输出的第一主成分包含最多的数据变异,第二主成分次之,之后的主成分包含的变异程度依次递减。SPSS输出结果如下:

    本研究中共有25个变量,那总特征值(eigenvalues of variance)是25,即每个变量自身的特征值为1。

    Total栏提示的是各主成分对数据变异的解释程度。

    以第一主成分为例,其特征值为6.730,占总体变异的6.730/25×100 = 26.919% (% of Variance栏)。同理,第二主成分的特征值为3.342,占总体变异的13.369%,以此类推。

    一般来说,如果某一项主成分的特征值小于1,那么我们就认为该主成分对数据变异的解释程度比单个变量小,应该剔除。本研究结果如下:

    第五主成分的特征值为1.049,大于1;而第六主成分的特征值为0.951,小于1,即应该保留前五位的主成分,剔除剩余部分。

    结论

    本研究采用主成分分析,通过25项问题调查315位应聘者的工作能力。

    研究变量之间存在线性相关关系(每组变量之间的相关系数均大于0.3),数据结构合理(KMO检验系数为0.833,单个变量的KMO检验系数均大于0.7,Bartlett's检验结果为P<0.001),提示研究数据可以进行主成分提取。< span="">

    主成分提取结果:研究提取前四位主成分。提取后的主成分累计解释59.9%的数据变异,分别反映应聘者的工作积极性、工作自主性、工作热情和工作责任感(如下图)

    展开全文
  • 是最外层的用户明细数据 数据清洗层主要是数据原始层的数据经过简单数据清洗之后的数据层, 主要是去除明显是脏数据, 比如年龄大于200岁, 地域来自 FFFF的 等明显异常数据 数据汇总层的数据主要是根据数据分析的...

    01

        写在前面

        

        

    我们经常在淘宝上购物, 作为淘宝方, 他们肯定想知道他的使用用户是什么样的, 是什么样的年龄性别, 城市, 收入, 他的购物品牌偏好, 购物类型, 平时的活跃程度是什么样的, 这样的一个用户描述就是用户画像分析

    无论是产品策划还是产品运营, 前者是如何去策划一个好的功能, 去获得用户最大的可见的价值以及隐形的价值, 必须的价值以及增值的价值, 那么了解用户, 去做用户画像分析, 会成为数据分析去帮助产品做做更好的产品设计重要的一个环节。

    那么作为产品运营, 比如要针用户的拉新, 挽留, 付费, 裂变等等的运营, 用户画像分析可以帮助产品运营去找到他们的潜在的用户, 从而用各种运营的手段去触达。

    因为当我们知道我们的群体的是什么样的一群人的时候, 潜在的用户也是这样的类似的一群人, 这样才可以做最精准的拉新, 提高我们的ROI

    在真正的工作中, 用户画像分析是一个重要的数据分析手段去帮助产品功能迭代, 帮助产品运营做用户增长。

    总的来说, 用户画像分析就是基于大量的数据,  建立用户的属性标签体系, 同时利用这种属性标签体系去描述用户

    02

       用户画像的作用

    像上面描述的那样, 用户画像的作用主要有以下几个方面

    1.广告投放

    在做用户增长的例子中, 我们需要在外部的一些渠道上进行广告投放, 对可能的潜在用户进行拉新, 比如B站在抖音上投广告

    我们在选择平台进行投放的时候, 有了用户画像分析, 我们就可以精准的进行广告投放, 比如抖音的用户群体是18-24岁的群体, 那么广告投放的时候就可以针对这部分用户群体进行投放, 提高投放的ROI

    假如我们没有画像分析, 那么可能会出现投了很多次广告, 结果没有人点击

    2.精准营销

    假如某个电商平台需要做个活动给不同的层次的用户发放不同的券, 那么我们就要利用用户画像对用户进行划分, 比如划分成不同的付费的活跃度的用户, 然后根据不同的活跃度的用户发放不用的优惠券。

    比如针对付费次数在 [1-10] 的情况下发 10 元优惠券刺激, 依次类推

    3. 个性化推荐

    精确的内容分发, 比如我们在音乐app 上看到的每日推荐, 网易云之所以推荐这么准, 就是他们在做点击率预估模型(预测给你推荐的歌曲你会不会点击)的时候, 考虑了你的用户画像属性。

    比如根据你是90后, 喜欢伤感的, 又喜欢杰伦, 就会推荐类似的歌曲给你, 这些就是基于用户画像推荐

    4. 风控检测

    这个主要是金融或者银行业设计的比较多, 因为经常遇到的一个问题就是银行怎么决定要不要给一个申请贷款的人给他去放贷

    经常的解决方法就是搭建一个风控预测模型, 去预约这个人是否会不还贷款,同样的, 模型的背后很依赖用户画像。

    用户的收入水平, 教育水平, 职业, 是否有家庭, 是否有房子, 以及过去的诚信记录, 这些的画像数据都是模型预测是否准确的重要数据

    5. 产品设计

    互联网的产品价值 离不开 用户 需求 场景 这三大元素, 所以我们在做产品设计的时候, 我们得知道我们的用户到底是怎么样的一群人, 他们的具体情况是什么, 他们有什么特别的需求, 这样我们才可以设计出对应解决他们需求痛点的产品功能

    在产品功能迭代的时候, 我们需要分析用户画像行为数据, 去发现用户的操作流失情况, 最典型的一种场景就是漏斗转化情况, 就是基于用户的行为数据去发现流失严重的页面, 从而相对应的去优化对应的页面,

    比如我们发现从下载到点击付款转化率特别低,那么有可能就是我们付款的按钮的做的有问题, 就可以针对性的优化按钮的位置等等

    同时也可以分析这部分转化率主要是在那部分用户群体中低, 假如发现高龄的用户的转化率要比中青年的转化率低很多, 那有可能是因为我们字体的设置以及按钮本身位置不显眼等等, 还有操作起来不方便等等因素

    6. 数据分析

    在做描述性的数据分析的时候, 经常需要画像的数据, 比如描述抖音的美食博主是怎么样的一群人, 他们的观看的情况, 他们的关注其他博主的情况等等

    简单来说就是去做用户刻画的时候, 用户画像可以帮助数据分析刻画用户更加清晰。

    03

    如何搭建用户画像

    用户画像搭建的架构如下: 

    数据层: 

    首先 是数据层,  用户画像的基础是首先要去获取完整的数据, 互联网的数据主要是 利用打点, 也就是大家说的数据埋点上报上来的, 整个过程就是 数据分析师会根据业务需要提数据上报的需求,然后由开发完成, 这样就有了上报的数据。

    除了上报的数据, 还有其他数据库同步的数据, 一般会把数据库的数据同步到hive表中, 按照数据仓库的规范, 按照一个个主题来放置

    还有一些其他的数据比如外部的一些调研的数据, 以excel 格式存在, 就需要把excel 数据导入到hive 表中

    挖掘层:

    有了基础的数据以后, 就进入到挖掘层, 这个层次主要是两件事情, 一个是数据仓库的构建, 一个是标签的预测, 前者是后者的基础。

    一般来说我们会根据数据层的数据表, 对这些数据表的数据进行数据清洗,数据计算汇总, 然后按照数据仓库的分层思想, 比如按照 数据原始层, 数据清洗层, 数据汇总层, 数据应用层等等进行表的设计

    数据原始层的表的数据就是上报上来的数据入库的数据, 这一层的数据没有经过数据清洗处理, 是最外层的用户明细数据

    数据清洗层主要是数据原始层的数据经过简单数据清洗之后的数据层, 主要是去除明显是脏数据, 比如年龄大于200岁,  地域来自 FFFF的 等明显异常数据

    数据汇总层的数据主要是根据数据分析的需求, 针对想要的业务指标, 比如用户一天的听歌时长, 听歌歌曲数, 听的歌手数目等等, 就可以按照用户的维度, 把他的行为进行聚合, 得到用户的轻量指标的聚合的表。

    这个层的用处主要是可以快速求出比如一天的听歌总数, 听歌总时长, 听歌时长高于1小时的用户数, 收藏歌曲数高于100 的用户数是多少等等的计算就可以从这个层的表出来

    数据应用层主要是面向业务方的需求进行加工, 可能是在数据汇总的基础上加工成对应的报表的指标需求, 比如每天听歌的人数, 次数, 时长, 搜索的人数, 次数, 歌曲数等等

    按照规范的数据仓库把表都设计完成后, 我们就得到一部分的用户的年龄性别地域的基础属性的数据以及用户观看 付费 活跃等等行为的数据

    但是有一些用户的数据是拿不到的比如音乐app 为例, 我们一般是拿不到用户的听歌偏好这个属性的数据, 我们就要通过机器学习的模型对用户的偏好进行预测

    机器学习的模型预测都是基于前面我们构建的数据仓库的数据的, 因为只有完整的数据仓库的数据, 是模型特征构建的基础

    服务层:

    有了数据层和挖掘层以后, 我们基本对用户画像体系构建的差不多, 那么就到了用户画像赋能的阶段。

    最基础的应用就是利用用户画像宽表的数据, 对用户的行为进行洞察归因 挖掘行为和属性特征上的规律

    另外比较大型的应用就是搭建用户画像的平台, 背后就是用户画像表的集成。

    用户提取: 我们可以利用用户画像平台, 进行快速的用户选取,  比如抽取18-24岁的女性群体 听过杰伦歌曲的用户, 我们就可以快速的抽取。

    分群对比: 我们可以利用画像平台进行分群对比。比如我们想要比较音乐vip 的用户和非vip 的用户他们在行为活跃和年龄性别地域 注册时间, 听歌偏好上的差异, 我们就可以利用这个平台来完成

    功能画像分析: 我们还可以利用用户画像平台进行快速进行某个功能的用户画像描述分析, 比如音乐app 的每日推荐功能, 我们想要知道使用每日推荐的用户是怎么样的用户群体, 以及使用每日推荐不同时长的用户他们的用户特征分别都是怎么样的,就可以快速的进行分析

    详解用户流失原因分析该如何入手?

    12000+字超详细 SQL 语法速成!

    后台回复“入群”即可加入小z干货交流群
    
    展开全文
  • 是一种旨在寻找隐藏在多变量数据中、无法直接观察到却影响或支配可测变量的潜在因子、并估计潜在因子对可测变量的影响程度以及潜在因子之间的相关性的一种多元统计分析方法 基本思想 根据相关性大小把变量...

     

    因子分析法是指从研究指标相关矩阵内部的依赖关系出发,把一些信息重叠、具有错综复杂关系的变量归结为少数几个不相关的综合因子的一种多元统计分析方法。 

    是一种旨在寻找隐藏在多变量数据中、无法直接观察到却影响或支配可测变量的潜在因子、并估计潜在因子对可测变量的影响程度以及潜在因子之间的相关性的一种多元统计分析方法

    基本思想

    根据相关性大小把变量分组,使得同组内的变量之间相关性较高,但不同组的变量不相关或相关性较低,每组变量代表一个基本结构一即公共因子。

    为什么做因子分析

    举例说明:在实际门店问题中,往往我们会选择潜力最大的门店作为领航店,以此为样板,实现业绩和利润的突破及未来新店的标杆。选择领航店过程中我们要注重很多因素,比如:

    ↘所在小区的房价

    ↘总面积

    ↘户主年龄分布

    ↘小区户数

    ↘门店面积

    ↘2公里范围内竞争门店数量等

    收集到所有的这些数据虽然能够全面、精准的确定领航店的入选标准,但实际建模时这些变量未必能够发挥出预期的作用。主要体现两方面:计算量的问题;变量间的相关性问题。

    这时,最简单直接的方案就是削减变量个数,确定主要变量,因子分析以最少的信息丢失为前提,将众多的原有变量综合成少数的综合指标。

    因子分析特点

    因子个数远小于变量个数;

    能够反应原变量的绝大数信息;

    因子之间的线性关系不显著;

    因子具有命名解释性

    因子分析步骤

    1.原有变量是否能够进行因子分析;

    2.提取因子;

    3.因子的命名解释;

    4.计算因子得分;五、综合评价

    因子与主成分分析的区别

    相同:都能够起到处理多个原始变量内在结构关系的作用

    不同:主成分分析重在综合原始变适的信息.而因子分析重在解释原始变量间的关系,是比主成分分析更深入的一种多元统计方法

    因子分析可以看做是优化后的主成分分析,两种方法有很多共通的地方,但应用方面各有侧重。

    因子分析应用场景

    因子分析方法主要用于三种场景,分别是:

    l信息浓缩:将多个分析项浓缩成几个关键概括性指标。比如将多个问卷题浓缩成几个指标。如果偏重信息浓缩且关注指标与分析项对应关系,使用因子分析更为适合。

    l权重计算:利用方差解释率值计算各概括性指标的权重。在信息浓缩的基础上,可进一步计算每个主成分/因子的权重,构建指标权重体系。

    l综合竞争力:利用成分得分和方差解释率这两项指标,计算得到综合得分,用于综合竞争力对比(综合得分值越高意味着竞争力越强)。此类应用常见于经济、管理类研究,比如上市公司的竞争实力对比。

    因子分析案例

    现在有 12 个地区的 5 个经济指标调查数据(总人口、学校校龄、总雇员、专业服务、中等房价),为对这 12 个地区进行综合评价,请确定出这 12 个地区的综合评价指标。(综合竞争力应用场景


    同一指标在不同地区是不同的,用单一某一个指标难以对12个地区进行准确的评价,单一指标只能反映地区的某一方面。所以,有必要确定综合评价指标,便于对比。因子分析方法就可以应用在这个案例中。

    5 个指标即为我们分析的对象,我们希望从这5个可观测指标中寻找出潜在的因素,用这些具有综合信息的因素对各地区进行评价。

    下图spss因子分析的操作界面主要包括5方面的选项,变量区只能选择数值型变量,分类型变量不能进入该模型。

    spss软件为了消除不同变量间量纲和数量级对结果的影响,在该过程中默认自动进行标准化处理,因此不需要对这些变量提前进行标准化处理。

    描述统计选项卡

    希望看到各变量的描述统计信息,要对比因子提取前后的方差变化,选定“单变量描述性”和“原始分析结果”;

    现在是基于相关矩阵提取因子,所以,选定相关矩阵的“系数和显著性水平“,

    另外,比较重要的还有 KMO 和球形检验,通过KMO值,我们可以初步判断该数据集是否适合采用因子分析方法,kmo结果有时并不会出现,这主要与变量个数和样本量大小有关。

    抽取选项卡:在该选项卡中设置如何提取因子

    提取因子的方法有很多,最常用的就是主成分法。

    因为参与分析的变量测度单位不同,所以选择“相关矩阵”,如果参与分析的变量测度单位相同,则考虑选用协方差矩阵。

    经常用到碎石图对于判断因子的个数很有帮助,一般都会选择该项。关于特征值,一般spss默认只提取特征值大于1的因子。收敛次数比较重要,可以从首次结果反馈的信息进行调整。

    因子旋转选项卡

    因子分析要求对因子给予命名和解释,是否对因子旋转取决于因子的解释。

    旋转就是坐标变换,使得因子系数向1 和 0 靠近,对公因子的命名和解释更加容易。旋转方法一般采用”最大方差法“即可,输出旋转后的因子矩阵和载荷图,对于结果的解释非常有帮助。

    如果不经旋转因子已经很好解释,那么没有必要旋转,否则,应该旋转。

    保存因子得分

    要计算因子得分就要先写出因子的表达式。因子是不能直接观察到的,是潜在的。但是可以通过可观测到的变量获得。

    因子分析模型是原始变量为因子的线性组合,现在我们可以根据回归的方法将模型倒过来,用原始变量也就是参与分析的变量来表示因子。从而得到因子得分。因子得分作为变量保存,对于以后深入分析很有用处。

    结果解读:验证数据是否适合做因子分析

    参考kmo结果,一般认为大于0.5,即可接受。同时还可以参考相关系数,一般认为分析变量的相关系数多数大于 0.3,则适合做因子分析;

    KMO=0.575 检验来看,不是特别适合因子分析,基本可以通过。

    结果解读:因子方差表

    提取因子后因子方差的值均很高,表明提取的因子能很好的描述这 5 个指标。

    方差分解表表明,默认提取的前两个因子能够解释 5 个指标的 93.4%。碎石图表明,从第三个因子开始,特征值差异很小。综上,提取前两个因子。

    结果解读:因子矩阵

    旋转因子矩阵可以看出,经旋转后,因子便于命名和解释。

    因子 1主要解释的是中等房价、专业服务项目、中等校平均校龄,可以命名为社会福利因子;

    因子 2 主要解释的是其余两个指标,总人口和总雇员。可以命名为人口因子。

    因子分析要求最后得到的因子之间相互独立,没有相关性,而因子转换矩阵显示,两个因子相关性较低。可见,对因子进行旋转是完全有必要的。

    结果解读:因子系数

    因子得分就是根据这个系数和标准化后的分析变量得到的。在数据视图中可以看到因子得分变量。

    结论

    经过因子分析实现了目的,找到了两个综合评价指标,人口因子和福利因子。

    从原来的 5 个指标挖掘出 2 个潜在的综合因子。可以对12 个地区给出客观评价。

    可以根据因子1或因子2得分,对这12个地区进行从大到小排序,得分高者被认为在这个维度上有较好表现。

    展开全文
  • 数据分析的基本思路 (1)从ncbi的geo或者其它数据库中查找自己感兴趣的RNASeq数据,至少要求给出如下信息: 该套数据所发表的文章的名字: 该套数据的下载网址: 该套数据基本情况介绍(简介以及该套数据包含...

    数据分析的基本思路

    (1)从ncbi的geo或者其它数据库中查找自己感兴趣的RNASeq数据,至少要求给出如下信息:

    (2)对芯片数据进行质量控制评价及处理(如果质量差的话,每个样本都应该处理),

      可以用软件Fastqc+Trimmomatic配合使用,也可以用其它软件替换

    (3)用TopHat2 + Cufflinks+Hisat系列软件进行分析

    (4)用Hisat、StringTie和Ballgown系列软件进行分析

    (5)(3)和(4)的结果都要进行差异表达分析和聚类分析

    (6)比较(3)和(4)两种系列软件结果的差异表达分析和聚类分析

    (一)数据集选择

    1、根据筛选标准(是否为RNA-seq数据、样本量、样本分类清晰程度、是否有sra原始数据),到GEO数据库中查找数据集。

    数据ID:GSE111845

    数据来源:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111845

    样本介绍:Standard methods were used for RNA-sequencing library construction, EXBead preparation, and Next-Generation sequencing, based on the protocol provided for the Life Technologies SOLiD4 system. Briefly, 2μg of total RNA per sample of 30 human liver samples (fetal, pediatric, and adult; n=10) were used for library preparation. The rRNA was depleted using RiboMinus Eukaryote Kit for RNA-Seq.

    分类:成人、幼儿、胎儿的肝样本各十个。

    sra accession number:SRP135703

    一共下载了10个样本,其中成人、胎儿各5个。

    SRR6835938

    幼儿

    SRR6835935

    成人

    SRR6835936

    成人

    SRR6835940

    幼儿

    SRR6835941

    幼儿

    SRR6835942

    幼儿

    SRR6835943

    幼儿

    SRR6835944

    成人

    SRR6835945

    成人

    SRR6835947

    成人

    (二)数据下载

    【参考链接】https://www.jianshu.com/p/8dca09077df3

    1、使用SRAtoolkitprefetch命令下载SRA数据。

    for i in `seq 35 47`; 
    do 
    prefetch SRR68359${i} 
    done
    

    2、使用SRAtoolkitfastq-dumpSRA数据解析成fasta文件。

    由于sra文件是单端测序(此处具体参考sratoolkit工具的操作说明),故而使用fastq-dump指令为:

    fastq-dump  SRR6835947.sra

    (三)质量控制和预处理

    1、使用fastqc(软件的安装步骤自行百度)对读段数据进行质控分析。

    fastqc  -f  fastq  -o  result  SRR6835935.out.fastq 

    fastqc的使用说明:

    fastqc是用来评估测序reads质量的软件。

    使用命令:fastqc -f fastq -o result reads.fq.gz

    #-f 表示的是 输入文件的类型

    #-o 表示的是 输出的目录(事先用mkdir创建了一个文件夹用于存放结果)

    #reads.fq 表示的是被评估的测序数据

    #在Linux平台敲入指令:
    s45@HP45:~/Biosofts/FastQC$ fastqc -f fastq -o result SRR6835935.out.fastq 
    
    #屏幕回显:
    Started analysis of SRR6835935.out.fastq
    Approx 5% complete for SRR6835935.out.fastq
    Approx 10% complete for SRR6835935.out.fastq
    Approx 15% complete for SRR6835935.out.fastq
    Approx 20% complete for SRR6835935.out.fastq
    Approx 25% complete for SRR6835935.out.fastq
    Approx 30% complete for SRR6835935.out.fastq
    Approx 35% complete for SRR6835935.out.fastq
    Approx 40% complete for SRR6835935.out.fastq
    Approx 45% complete for SRR6835935.out.fastq
    Approx 50% complete for SRR6835935.out.fastq
    Approx 55% complete for SRR6835935.out.fastq
    Approx 60% complete for SRR6835935.out.fastq
    Approx 65% complete for SRR6835935.out.fastq
    Approx 70% complete for SRR6835935.out.fastq
    Approx 75% complete for SRR6835935.out.fastq
    Approx 80% complete for SRR6835935.out.fastq
    Approx 85% complete for SRR6835935.out.fastq
    Approx 90% complete for SRR6835935.out.fastq
    Approx 95% complete for SRR6835935.out.fastq
    Analysis complete for SRR6835935.out.fastq
    

     

    2、质控分析结果

     

    A

    B

     

     

    C

    D

     

     

    E

    F

     

     

    说明:这篇博客主要是介绍RNA-seq总的分析流程。关于fastqc所得到的这许多结果图的生物学意义,稍后会另外写一篇文章详细介绍。

    3、使用Trimmomatic软件对读段数据进行过滤

    java -jar trimmomatic-0.38.jar SE -phred33 ./RNA-seq-data/$ip ./OUTresult/ SRR6835935.out.fastq.gz ILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 HEADCROP:10 SLIDINGWINDOW:4:15 MINLEN:36

    注明:本次操作我所使用的序列文件,在fastqc评估时评估的各项指标都较好,认为此处不需要对数据再进行过滤,使原始数据失真。各位在实际的操作的过程,应该根据上面那一步的fastqc的结果来进行判断是否trim,以及trim哪里,trim多少。

     

    (四)将读段比对到参考基因组

    【准备工作】

    1)人类基因组序列文件的下载

    下载NCBI的genome数据库中下载人类基因组文件(GRCh37/hg19,根据自己的实际情况进行选择,选择错误的话,下游的操作会报错)

    ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh37_latest/refseq_identifiers/GRCh37_latest_genomic.fna.gz

    2)人类基因组注释文件的下载

    在GENECODE数据库中下载chr开头的gff3人类基因组注释文件。

    https://www.gencodegenes.org/human/release_29.html

    本次实验下载hg19类型的Comprehensive gene annotation(Regions:CHR)人类染色体的gff3格式的注释文件。

     

    【方案一】

    1、使用TopHat程序将读段比对到参考基因组

    (1)构建参考索引

    nohup bowtie2-build -f  human_genomic.fna human &

    #nohup ……&是将程序挂在服务器后台运行的指令。具体可自行百度。

    2)将读段比对到基因组

     

    for i in *fastq.gz
      do
      echo $i
         tophat2 -p 4 -o /home/s45/index/result/  /home/s45/index/human19   /home/s45/index/data/$i
      done
    

    #同样地,循环指令。根据自己的实际情况,修改参数。

    # 本次操作输出文件

    1.accepted_hits.bam 包含比对,为bam格式,比对根据染色体坐标被排序。

    2.junctions.bed 包含发现外显子结界,为bed格式。结界由两块组成,其中每个块与跨越该结界的任何读段的最大延伸量一样长。得分数是跨越结界的比对的数目。

    3.insertions.bed 包含发现的插入。

    4.deletion.bed 包含发现的缺失。

    5.align_summary.txt 报告

    【方案二】

    2、使用HiSat程序将读段比对到参考基因组

    (1)构建参考索引

    hisat-build  human19.fa  human3

    (2)将读段比对到基因组

    for i in `seq 41 47`
       do
             hisat2 -t --dta -p 8 -x human3 -U ./data/SRR68359${i}.out.fastq.gz  -S ./result/SRR68359${i}.sam
       done
    

    (五)转录组组装

    【方案一】

    1、使用cufflinks系列软件对读段进行组组装

    (1)cufflinks组装

    cufflinks -p 8 -o ./   35_accepted_hits.bam 

    参数说明:

    #-p 指定线程数

    #-o :指定输出文件所在目录

    #最后是Tophat比对后的结果文档bam文件

    2cuffmerge转录本合并

    cuffmerge -g genes.gtf  –s  ./genome.fa  -p 8   assemblies.txt

    # assemblies.txt文件的具体内容(指明要比对的gtf文件的路径)

    ./35/ transcripts.gtf
    ./36/ transcripts.gtf
    ./38/ transcripts.gtf
    ./40/ transcripts.gtf
    ./41/ transcripts.gtf
    ./42/ transcripts.gtf
    ./43/ transcripts.gtf
    ./44/ transcripts.gtf
    ./45/ transcripts.gtf
    ./47/ transcripts.gtf

    3cuffquant基因和转录本表达定量

    cuffquant -o  sample_quant  -p  8  35_accepted_hits.bam

    这些步骤生成的文件:

    【方案二】

    2、使用stringtie软件对读段进行组装

    (1)使用samtools将HiSeq产生的.sam文件转换为.bam文件

    (2)使用samtools将.bam文件进行排序

    samtools sort SRR6835936.bam SRR6835936.sort

    (3)使用stringtie对读段进行组装

    stringtie -p 8 -G ./hg19_refSeq.gtf -o ./gtf/36.gtf - 36 SRR6835936.sort.bam
    stringtie merge -p 8 -G ./genes/chr.gtf -o stringtie_merged.gtf ./mergelist.txt

    # mergelist.txt的具体内容(类似于前面的某一步):

    ./35.gtf
    ./36.gtf
    ./38.gtf
    ./40.gtf
    ./41.gtf
    ./42.gtf
    ./43.gtf
    ./44.gtf
    ./45.gtf
    ./47.gtf
    
    stringtie -e -B -p 8 -G ./gtf/merged.gtf -o ./ballgown/35/365gtf  SRR6835935.sort.bam

    #e2t.ctab    外显子和转录本间的关系

    #e_data.ctab

    #i2t.ctab    内含子和转录本间的关系

    #i_data.ctab

    #t_data.ctab  转录本的数据

    (六)差异表达分析

    【方案一】

    1、使用Hiseq软件对TopHat比对结果进行定量

    htseq-count -f bam /media/zxx/6476-C43F2/SRR6835936_accepted_hits.bam /media/zxx/6476-C43F2/hg19_refSeq.gtf >383_sample.count

    2、使用R语言,将Hiseq定量结果,汇总成表达矩阵

    setwd("D:/Rworkplace/HTSeq")
    Adult1 <- read.table("D:/Rworkplace/HTSeq/35_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Adult2 <- read.table("D:/Rworkplace/HTSeq/36_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Pediatric1<-read.table("D:/Rworkplace/HTSeq/38_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Pediatric2<-read.table("D:/Rworkplace/HTSeq/40_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Pediatric3<-read.table("D:/Rworkplace/HTSeq/41_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Pediatric4<-read.table("D:/Rworkplace/HTSeq/42_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Pediatric5<-read.table("D:/Rworkplace/HTSeq/43_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Adult3<-read.table("D:/Rworkplace/HTSeq/44_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Adult4<-read.table("D:/Rworkplace/HTSeq/45_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    Adult5<-read.table("D:/Rworkplace/HTSeq/47_sample.count", sep="\t", col.names = c("gene_id","readsSum"))
    
    myMatrix <- merge(merge(Adult1, Adult2,Adult3,Adult4,Adult5,by="gene_id"), merge(Pediatric1,Pediatric2,Pediatric3,Pediatric4,Pediatric5,by="gene_id"))
    
    Merge_func<-function(x,y){
      df <- merge(x,y,by="gene_id")
      rownames(df)<- df$Row.names
      df$Row.names<-NULL
      return(df)
    }
    Adult<- Reduce(Merge_func,list(Adult1, Adult2,Adult3,Adult4,Adult5))
    Pediatric<-Reduce(Merge_func,list(Pediatric1,Pediatric2,Pediatric3,Pediatric4,Pediatric5))
    myMatrix<-merge(Adult,Pediatric,by="gene_id")
    colnames(myMatrix)<-c("gene_id","Adult_35","Adult_36","Adult_44","Adult_45","Adult_47","Pediatric_38","Pediatric_40","Pediatric_41","Pediatric_42","Pediatric_43")
    ENSEMBL <- gsub("\\.\\d*", "",myMatrix$gene_id)
    head(ENSEMBL, n=7)
    myMatrixS<-cbind(ENSEMBL,myMatrix)
    row.names(myMatrixS)<-myMatrixS[,1]
    myMatrixR<-myMatrixS[,-c(1:2)]
    myMatrixE<-myMatrixR[-c(1:5),]
    
    #myMatrixE就是最后的表达谱矩阵
    write.table(myMatrixE,file="HTSeqmatrix.txt",col.names = T,row.names = T)
    
    

    构建出的表达矩阵如下:

    "Adult_35" "Adult_36" "Adult_44" "Adult_45" "Adult_47" "Pediatric_38" "Pediatric_40" "Pediatric_41" "Pediatric_42" "Pediatric_43"

    "ENSG00000000003" 5 0 0 0 0 0 3 0 1 2

    "ENSG00000000005" 0 0 0 1 0 0 0 0 0 0

    "ENSG00000000419" 0 0 0 0 0 0 0 0 2 0

    "ENSG00000000457" 38 10 36 26 5 19 32 24 51 44

    "ENSG00000000460" 33 3 9 10 0 6 11 7 7 13

    "ENSG00000000938" 0 0 0 0 0 1 0 0 0 0

    "ENSG00000000971" 626 216 289 209 140 202 270 0 606 345

    "ENSG00000001036" 1 0 0 0 0 0 0 0 0 0

    "ENSG00000001084" 5 1 0 3 0 0 0 0 3 2

    "ENSG00000001167" 12 2 5 5 0 1 0 0 3 3

    "ENSG00000001460" 0 0 0 0 0 0 0 0 0 0

    "ENSG00000001461" 0 0 0 0 0 0 0 0 1 0

    "ENSG00000001497" 0 0 0 1 0 0 0 0 0 0

    "ENSG00000001561" 0 0 0 0 0 0 0 0 0 0

    "ENSG00000001617" 21 4 28 19 5 29 45 0 25 23

    "ENSG00000001626" 0 0 0 0 0 0 0 0 2 0

    "ENSG00000001629" 0 0 0 0 0 0 0 0 0 0

    "ENSG00000001630" 7 5 2 3 0 5 1 0 4 0

    3、对表达矩阵进行差异表达分析,筛选出差异基因

    Pvalue<0.05&log2(Foldchange)>1为阈值,筛选差异表达基因。

    setwd("D:/Rworkplace/HTseq ")
    geneexprSet<-read.table("microRNA.txt",header = T,sep = " ")
    dim(geneexprSet)
    TTest<-function(geneexprSet,numStart,numEnd){
      P_value<-c()
      Foldchange<-c()
      len<-dim(geneexprSet)[1]
      for(n in 1:len){
        controlsample<-as.numeric(geneexprSet[n,2:numStart])
        testsample<-as.numeric(geneexprSet[n,(numStart+1):numEnd])
        res<-t.test(controlsample,testsample,paired=TRUE)
        meanControl<-mean(controlsample)
        meanSample<-mean(testsample)
        fold<-meanSample/meanControl
        p<-res$p.value
        P_value<-append(P_value,p)
        Foldchange<-append(Foldchange,fold)
      }
      finalexprSet<-cbind(geneexprSet,Pvalue=P_value,Foldchange=Foldchange)
      write.table(finalexprSet,file="microfinalexprSet.txt",col.names = T,row.names = T)
    }
    TTest(geneexprSet,11,21)
    geneexprSet<-read.table("exprSet.txt",header = T,sep = " ")
    geneexprSet<-geneexprSet[order(geneexprSet$Pvalue), ]
    diffresult<-subset(geneexprSet,geneexprSet$Pvalue<0.05&log2(geneexprSet$Foldchange)>1,select = c(-Pvalue,-Foldchange))
    write.table(diffresult,file="diffresult.txt",col.names = T,row.names = T)
    

    筛选得到的差异表达矩阵如下所示:

    "Adult_35" "Adult_36" "Adult_44" "Adult_45" "Adult_47" "Pediatric_38" "Pediatric_40" "Pediatric_41" "Pediatric_42" "Pediatric_43"

    "ENSG00000023839" 0 1 0 0 0 2 3 0 1 3

    "ENSG00000031698" 18 14 19 22 7 37 23 19 42 47

    "ENSG00000048462" 0 0 0 0 0 1 1 0 1 1

    "ENSG00000066294" 4 1 4 1 0 2 9 7 7 8

    "ENSG00000085999" 0 0 0 0 0 4 2 1 3 4

    "ENSG00000086015" 394 236 304 552 276 522 576 610 849 1139

    4、对差异表达基因可视化

    (1)对矩阵进行预处理

    dealdata<- read.table("diffresult.txt",header = TRUE)
    #去除NA
    dealdataNA<-na.omit(dealdata)
    #在dealdata矩阵的基础上,绘图
    

     

    (2)绘制火山图查看全部基因的表达情况

    library(ggplot2)
    volcano<-subset(dealdataNA,select = c(Pvalue,Foldchange))
    threshold<-as.factor((log2(volcano$Foldchange)>1.5|log2(volcano$Foldchange)<(-1.5))&volcano$Pvalue<0.05)
    r03=ggplot(volcano,aes(log2(Foldchange),-log2(Pvalue),colour=threshold))+geom_point()+xlim(-10,10)
    r04=r03+geom_vline(xintercept=c(-1.5,1.5),linetype="dotted",size=1)+geom_hline(yintercept=-log2(0.05),col="blue")
    r05=r04+labs(title="Volcanoplot")+theme(plot.title = element_text(hjust = 0.5))
    

    (3)对筛选出的差异表达基因绘制热图

    diffresult<-subset(dealdataNA,dealdataNA$Pvalue<0.05&log2(dealdataNA$Foldchange)>1,select = c(-Pvalue,-Foldchange))
    write.table(diffresult,file="diffresult.txt",col.names = T,row.names = T)
    library(pheatmap)
    bk = unique(c(seq(-1,1, length=100)))
    col_anno<-data.frame(type=c(rep("adult",5),rep("pediatric",5)),row.names=colnames(diffresult))
    
    pheatmap(diffresult,
             breaks=bk,
             main="the heatmap for diffgene expression in age samples(P<0.05)",
             scale='row',
             cutree_rows=2,
             cutree_cols=4,
             annotation_col=col_anno,
             fontsize_col=15,
             treeheight_row=120,
             treeheight_col=20,
             show_rownames = F ,
             color = colorRampPalette(c("red","black","green"))(100))
    

    对差异表达基因进行计数,一共筛选出来85个差异表达基因(不分上调/下调)。

    【方案二】

    1、安装R包Ballgown包

    #Ballgoon包的安装
    BiocManager::install("Ballgoon", ask=F, update=F)
    #安装成功!
    #其余genefiliter/dplyr/devtools通过常规方法可以完成。
    #其次安装RSkittleBrewer
    library(devtools)
    install_github('alyssafrazee/RSkittleBrewer')
    #安装成功!
    #加载包
    library(RSkittleBrewer)
    library(ballgown)
    library(dplyr)
    library(genefilter)
    library(devtools)
    

    2、使用R语言的程序包Ballgown对Stringtie的组装结果提取分析

    setwd("D:\\Rworkplace\\ballgown")
    sample<-read.csv("sample.csv")
    bg = ballgown(dataDir='./', samplePattern='SRA', pData=sample)
    bg_filt=subset(bg,"rowVars(texpr(bg))>1",genomesubset=TRUE)
    result_transcripts=stattest(bg_filt,feature = "transcript",covariate="sample", getFC=TRUE,meas="FPKM")
    result_genes=stattest(bg_filt,feature = "gene",covariate = "sample", getFC=TRUE,meas="FPKM")
    result_transcripts=data.frame(geneNames=ballgown::geneNames(bg_filt),geneIDs=ballgown::geneIDs(bg_filt),result_transcripts)
    result_transcripts=arrange(result_transcripts,pval)
    result_genes=arrange(result_genes,pval)
    write.csv(result_transcripts, "chrX_transcript_results.csv",row.names=FALSE)
    write.csv(result_genes, "chrX_gene_results.csv",row.names=FALSE)
    

    3、差异表达分析

    difftranscript<-subset(result_transcripts,result_transcripts$qval<0.05)
    diffgene<-subset(result_genes,result_genes$qval<0.05)
    tropical= c("green","red")
    palette(tropical)
    fpkm = texpr(bg,meas="FPKM")
    fpkm = log2(fpkm+1)
    
    setwd("D:\\Rworkplace\\ballgown")
    gene<-read.csv("chrX_gene_results.csv")
    
    library(ggplot2)
    volcano<-subset(gene,select = c(pval,fc))
    threshold<-as.factor((log2(volcano$fc)>1.5|log2(volcano$fc)<(-1.5))&volcano$pval<0.05)
    tropical= c("green","red")
    palette(tropical)
    r03=ggplot(volcano,aes(log2(fc),-log2(pval),colour=threshold))+geom_point()+xlim(-10,10)
    r04=r03+geom_vline(xintercept=c(-1.5,1.5),linetype="dotted",size=1)+geom_hline(yintercept=-log2(0.05),col="blue")
    r05=r04+labs(title="Volcanoplot")+theme(plot.title = element_text(hjust = 0.5))
    

    4、可视化

     (1)绘制盒氏图

    boxplot(fpkm,col=as.numeric(sample$sample),las=2,ylab='log2(FPKM+1)',cex.axis=0.7,ylim=c(0,8),main="FPKM, adult in green, pediatric in red")

      

    这个图归一化了?我是不信的。时间原因,就不重新做了。大家在操作的过程中注意。

    (2)转录本分布图

    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835935'),samples="SRA35")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835936'),samples="SRA36")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835938'),samples="SRA38")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835940'),samples="SRA40")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835942'),samples="SRA42")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835943'),samples="SRA43")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835944'),samples="SRA44")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835945'),samples="SRA45")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835947'),samples="SRA47")
    plotTranscripts(ballgown::geneIDs(bg)[21808], bg, main=c('transcripts MSTRG.11191 in sample SRR6835941'),samples="SRA41")
    

    以差异转录本MSTRG.11191为例,观察各样本在转录本上的分布。

    Adult

    Pediatric

    A

    F

    B

    G

    C

    H

    D

    I

    E

    J

     

    plotMeans('MSTRG.11191', bg_filt,groupvar="sample",legend=T)
    不同的R包的出现,一定是为了解决特定的需求。这在我们使用的过程中是尤其要注意的。

     

    比较转录本MSTRG.11191在不同的样本组之间表达明显不同。

    plot(fpkm[21808,]~sample$sample, border=c(1,2), main=paste(ballgown::geneIDs(bg)[21808],' : ', ballgown::transcriptNames(bg)[21808]),pch=19, xlab="age",ylab='log2(FPKM+1)',cex=1)
    points(fpkm[21808,]~jitter(as.numeric(sample$sample)),col=c("black","blue"),pch=c(19,20))
    grid()
    

    (3)火山图

    筛选得到309个差异表达基因(不分上下调)。

    (七)分析与比较(相对简单)

     

    参考链接:https://blog.csdn.net/hill_night/article/details/44829965

    【HISAT与TopHat的对比】

    速度:100bp,20M的reads

    tophat用时26.7分钟,tophat2 1170分钟

    比对正确率:hisat(99.2%),tophat2(97.4%)

    敏感度和准确性:HISAT(97.3%,94.8%) Tophat2(90.6%,82.6%)

     

    【StringTie和Cufflinks对比】

    算法:

    cufflinks parsimony算法  (简约算法):生成最少的亚型

    说这种算法没有考虑转录丰度,在isoforms方面算的不准。其在算表达量的时候,按照图上的说法是用了最大似然冗余算法。

    stringTie先将reads分为不同的类,然后再针对每个类的reads生成一个拼接图来确定转录本,之后每个转录本产生一个流神经网络的最大流算法来评估表达水平

    这个算法的意思对应过来就是在一个基因处的若干个转录本,如何分配reads的数目才能让每个转录本的数目都处在最多的状态。

    组装结果:

    在组装方面StringTie具有一些优势,在低表达的部分,阈值过滤5%的StringTie比阈值过滤10%的准确度和敏感度还要高。

    关于组装效果,StringTie要好于cufflinks,同时他们又远远好于其他软件。

    找出来的基因中,cufflink找出来的70%在StringTie中有重合,相比于cufflink,StringTie在基因重构方面对三种类型的基因更有效,分别是:低冗余,高exon数目,和多重转录本。

    时间:

    对于100bp的reads而言

    StringTie 30min  cufflink 81min更多

     

    【cuffdiff和ballgown的对比】

    测试内存而言,cuffdiff  9个小时用148G,而ballgown 18秒用8G

    差异结果而言,cuffdiff更保守一些

    (八)总结与心得

    实验过程中遇到的错误及其解决方案

    1、TopHat运行组装程序时:

    Error:Could not find Bowtie 2 index files

    (/home/s45/Downloads/index/HomoGRCh38.96.4.*.bt2)

    (1)尝试在ensembel上下载同一来源的数据文件(gtf,fasta),建立索引。

    https://blog.csdn.net/sinat_38163598/article/details/73067785

    如果是在不同网站上下载的染色体序列和基因gtf注释文件,就可能会出现错误。

    所以下载数据时需要在同一网址下载,都在NCBI或者UCSC下载。

    尝试在ensembel上下载同一来源的数据文件(gtf,fasta),建立索引。

    依旧出错。

    tophat2-G/home/s45/Downloads/index/Homo_sapiens.GRCh38.96.4.gtf --transcriptome-index=GRCh38.96.4.tr GRCh38.96.4

    屏显:

    [2019-05-15 12:48:06] Building transcriptome files with TopHat v2.1.0
    -----------------------------------------------
    [2019-05-15 12:48:06] Checking for Bowtie
    		  Bowtie version:	 2.2.6.0
    [2019-05-15 12:48:16] Checking for Bowtie index files (genome)..
    Error: Could not find Bowtie 2 index files (GRCh38.96.4.*.bt2)

     

    (2)修改路径,依旧不行。

    tophat2  -G /home/s45/Downloads/index/Homo_sapiens.GRCh38.96.4.gtf --transcriptome-index=GRCh38.96.4.tr /home/s45/Downloads/index/GRCh38.96.4

     

    屏显:

    [2019-05-15 12:49:28] Building transcriptome files with TopHat v2.1.0
    -----------------------------------------------
    [2019-05-15 12:49:28] Checking for Bowtie
    		  Bowtie version:	 2.2.6.0
    [2019-05-15 12:49:36] Checking for Bowtie index files (genome)..
    Error: Could not find Bowtie 2 index files (/home/s45/Downloads/index/GRCh38.96.4.*.bt2)

     

    (3)尝试从tophat官网上index的下载:

    http://ccb.jhu.edu/software/tophat/igenomes.shtml

    依旧不行。

    最终解决方案:

    tophat2 -p 4 -o ./ /home/s45/index/human19  SRR6835935.fastq.gz

    2、Stringtie组装时出错。

    Error: the input alignment file is not sorted!

    我突然想到是不是需要建立索引。

    于是,回到samtools那一步,如何对文件进行排序。

    对sort后的bam文件建立索引。

    samtools index SRR6835935_sort_hisat.bam

    屏显:

    [E::hts_idx_push] NO_COOR reads not in a single block at the end 7 -1

     

    通过index指令,可以检验samtools文件是否经过排序。从上述的屏显,我可以推测出samtools对bam文件进行排序时,并没有排序成功。

    于是,重新进行排序。

    samtools sort SRR6835936.bam SRR6835936.sort

    屏显:

    s24@HP24:~/index/result$ samtools sort SRR6835936.bam SRR6835936.sort
    [bam_sort_core] merging from 6 files…

    再使用sort指令进行验证,这个sort文件能够建立索引,则推测这个文件是经过排序的。

    这种情况下,我们再次尝试代码。

    stringtie -p 8 -G ./hg19_refSeq.gtf -o ./36.gtf -l 36 SRR6835936.sort.bam

     

    在文件夹中发现36.gtf文档。

    推测可能原因是使用批处理指令时出错。

     

    3、使用Hiseq处理bam文件时报错。

    Error occured when reading beginning of SAM/BAM file.no BGZF EOF marker; file may be truncated

    对bam文件是否完整的诊断:

    samtools  view  42_align_sorted.bam|tail

    4、HiSeq定量结果全为零。

    ENSG00000000003    0
    ENSG00000000005    0
    ENSG00000000419    0
    ENSG00000000457    0
    ENSG00000000460    0
    ENSG00000000938    0
    ENSG00000000971    0
    ENSG00000001036    0
    ENSG00000001084    0
    ENSG00000001167    0
    ENSG00000001460    0
    ENSG00000001461    0
    ENSG00000001497    0
    ENSG00000001561    0
    ENSG00000001617    0
    ENSG00000001626    0
    ENSG00000001629    0
    ENSG00000001630    0
    ENSG00000001631    0
    ENSG00000002016    0

     

    gtf文件格式与tophat构建索引时的文件格式(1 or chr1)不一致。还有注意使用的注释文件的版本号。

    在比对、组装等操作时一定要注意注释文档、序列文件、索引文件的文件格式和版本的一致性。否则很有可能会出现问题。

    解决方案:

    在GENECODE数据库中可以下载到chr开头的hg19 版本的gff3格式的人类基因组注释文件。

    https://www.gencodegenes.org/human/release_29.html

    本次实验我主要下载的Comprehensive gene annotation(Regions:CHR)人类染色体的注释文件。

     

    【参考链接】(部分):

    https://www.cnblogs.com/chenpeng1024/p/9248891.html

    http://www.biotrainee.com/thread-415-1-1.html

    http://blog.sciencenet.cn/blog-759995-990471.html

    https://www.jianshu.com/p/681e02e7f9af

    https://blog.csdn.net/qq_40280759/article/details/80552193

    https://blog.csdn.net/weixin_34281477/article/details/87016243

    https://www.jianshu.com/p/1f5d13cc47f8?from=timeline

    https://www.jianshu.com/p/7ce1add65ab8

    http://blog.sina.com.cn/s/blog_9c93f6450102wcx9.html

    https://www.plob.org/article/12130.html

    http://blog.sciencenet.cn/blog-1509670-848266.html

    https://blog.csdn.net/kongxx/article/details/88653970

    http://www.360doc.com/content/18/0104/10/19913717_718924981.shtml

    https://www.jianshu.com/p/681e02e7f9af

    https://www.jianshu.com/p/796b5e711a26

    http://blog.sina.com.cn/s/blog_83f77c940102v7wl.html

    http://blog.sina.com.cn/s/blog_6836e9f40102v4hb.html

    https://www.cnblogs.com/wangprince2017/p/9937495.html

    https://blog.csdn.net/hill_night/article/details/44829965

    https://www.jianshu.com/p/38c2406367d5

    展开全文
  • 描述统计是通过图表或数学方法,对数据资料进行整理、分析,并对数据的分布状态、数字特征和随机变量之间关系进行估计和描述的方法。 描述统计分为集中趋势分析和离中趋势分析和相关分析三大部分。 △集中趋势分析 ...
  • 欢迎关注”生信修炼手册”!PCA是Principal components analysis的简称,叫做主成分分析,是使用最广泛的降维算法之一。所谓降维,就是降低特征的维度,最直观的变化就...
  • 1. 常用的多元统计分析方法有哪些? (1)多元正太分布检验 (2)多元方差-协方差分析 (3)聚类分析 (4)判别分析 (5)主成分分析 (6)对应分析 (7)因子分析 (8)典型相关性分析 (9)定性数据模型分析 (10...
  • 数据分析:常用的降维方法

    千次阅读 2021-04-23 14:46:14
    在统计学中,主成分分析是一种简化数据集的技术。它是一个线性变换。这个变换把数据变换到一个新的坐标系统中,使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上,第二大方差在第二个坐标(第二主成分)上...
  • MATLAB应用——数据分析与统计

    千次阅读 2021-04-22 02:18:14
    数学建模是用数学方法解决各种实际问题的桥梁,它已经渗透到各个领域,而且发挥出越来越重要的作用。面对自然科学和工程应用中的难题,大部分人无从入手,而个别人却能短时间内给出切实可行的解决方案,其差别往往...
  • Python 数据分析入门案例讲解 https://www.bilibili.com/video/BV18f4y1i7q9/ 1.标准化:去均值,方差规模化 Standardization标准化:将特征数据的分布调整成标准正态分布,也叫高斯分布,也就是使得数据的均值为...
  • 数据分析案例明确分析的目的数据处理原始数据数据清洗选择子子集重复数据处理缺失数据处理数据转化数据提取-字段分割异常值处理数据分析1.需求在哪里?2.需要什么样的人才?3.什么阶段需求最旺?结论 明确分析的目的...
  • 一种基于MATLAB的数据一致性的分析方法,其特征在于,包括以下步骤:1)准备好需要分析处理的原始数据;2)在MATLAB软件中利用GUIDE提供的空间设计生成GUI界面,并生成文件框架M;3)文件框架M将GUI界面初始化并设置回...
  • 16种常用的数据分析方法汇总* B2 M6 r( }8 X0 F" ~经常会有朋友问到一个朋友,数据分析常用的分析方法有哪些,我需要学习哪个等等之类的问题,今天数据分析精选给大家整理了十六种常用的数据分析方法,供大家参考...
  • 通过客户流失预测案例感悟数据分析设计方法,正如Gartner于2020年给出数据分析领域的技术趋势,更智能、更高速、更负责的AI,凸显新技术引领业务,数据驱动、AI驱动,以站在高纬度上的预测结果为顶层设计,倒逼数据...
  • 一、主成份分析模型图 二、主成份分析公式 三、主成份分析与因子分析步骤 ●数据预处理; ●选择分析模型; ●判断要选择的主成分/因子数目; ●选择主成分/因子; ●旋转主成分/因子; ●解释结果; ●计算主成分或因子...
  • 毕业设计 - 天气数据分析

    千次阅读 多人点赞 2021-10-08 18:40:15
    文章目录1 前言2 项目简介3 开始分析3.1 海洋对当地气候的影响3.2 导入数据集3.3 温度数据分析3.4 湿度数据分析3.5 风向频率玫瑰图3.6 计算风速均值的分布情况4 最后-毕设帮助 1 前言 Hi,大家好,这里是丹成学长,...
  • XPS原始数据处理(含分峰拟合)

    千次阅读 2021-01-14 16:23:09
    今天咱们更进一步,从测试中得到的原始数据开始谈起,分享一下XPS原始数据的处理,希望对大家有所帮助。前期内容:在前面两期分享中,我们都提到XPS定性分析包括全谱扫描和高分辨谱,其中全谱扫描不需要经过特别的...
  • 数据分析 数据规约

    千次阅读 2021-01-29 22:17:19
    一.概念 "数据规约"(Data Reduction)是指在尽可能保持数据原貌的前提下,最大限度地精简数据集.数据规约又分为2类:"属性规约"和"数值规约" 二....1.概念: ...降维方法有2大 类:"特征选择"和"特征提取
  • 数据分析 | 异常数据识别小结

    千次阅读 2021-03-16 15:31:53
    单变量数据异常识别2.1 简单统计量分析2.2 三倍标准差2.3 box-cox转化+3倍标准差基本介绍基本公式Box-Cox优势python 实现2.4 箱线图3. 时间序列数据异常识别3.1 设置恒定阈值3.2 设置动态阈值-移动平均法3.3 STL...
  • 5个步骤,用SPSS进行数据分析

    千次阅读 2020-12-28 21:30:16
    原标题:5个步骤,用SPSS进行数据分析 SPSS是一款非常强大的数据处理软件,那么该如何用SPSS进行数据分析呢?什么是SPSSSPSS是社会统计科学软件包的简称, 其官方全称为IBM SPSS Statistics。SPSS软件包最初由SPSS ...
  • 数据缺失值的4种处理方法

    千次阅读 2021-06-24 09:38:37
    机械原因是由于机械原因导致的数据收集或保存的失败造成的数据缺失,比如数据存储的失败,存储器损坏,机械故障导致某段时间数据未能收集(对于定时数据采集而言)。 人为原因是由于人的主观失误、历史局限或有意...
  • RNAseq数据分析--RNA-Seq数据质控

    千次阅读 2020-12-29 17:20:43
    数据分析之前,需要对数据质量控制 数据质控指标 碱基含量分布(应该满足碱基互补配对) 碱基质量分布 质量值>=Q20 : 好碱基 质量值<Q20: 坏碱基 测序质量软件 测序数据处理 adapter接头 ...
  • 数据分析实战

    千次阅读 2021-03-23 15:14:35
    数据分析实战数据分析基础数据分析全景图及修炼指南学习数据挖掘的最佳路径学数据分析要掌握哪些基本概念用户画像:标签化就是数据的抽象能力数据采集:如何自动化采集数据数据采集:如何用八爪鱼采集微博上的“D&...
  • XPS原始数据处理之 Avantage 软件篇

    千次阅读 2021-07-23 14:30:31
    其中全谱扫描不需要经过特别的处理,因此在这里不做过多介绍,大家对照前两期内容应该很容易实现全谱分析,得到自己想要的信息(主要是样品所含元素种类的鉴定),这里主要介绍高分辨谱的数据处理与分析。1XPS高分辨谱...
  • 认识工业大数据 什么是工业大数据?...工业大数据会配合工业互联网的技术,利用原始资料来支援管理上的决策。” 百度百科是这样说的:“工业大数据是指在工业领域中,围绕典型智能制造模式,从客户需
  • 常用SPSS数据处理方法,你都会吗?

    千次阅读 2020-12-28 21:30:16
    数据处理是在统计和分析数据时,第一步要做的。尤其是当面对大量数据时,数据处理是一个重要的过程,可以达到提高处理效率及精度的目的。为配合进行更好的分析,研究过程过可能涉及到以下数据处理工作:定义变量名、...
  • 数据分析中缺失值的处理方法1、缺失值的分类按照数据缺失机制可分为:(1)完全随机缺失(missing completely at random, MCAR)所缺失的数据发生的概率既与已观察到的数据无关,也与未观察到的数据无关.(2)随机缺失...
  • Hadoop大数据综合案例4-Hive数据分析

    千次阅读 2021-05-17 20:09:30
    数据分析是指用适当的统计分析方法对收集来的大量数据进行分析,从行业角度看,数据分析是基于某种行业目的,有目的的进行收集、整理、加工和分析数据的过程,通过提取有用信息,从而形成相关结论,这一过程也是质量...
  • 数据分析的过程描述

    千次阅读 2021-03-24 15:53:44
    数据分析过程可以用以下几步来描述: 转换和处理原始数据,以可视化方式呈现数据,建模做预测。 因此数据分析几乎可以概括为由以下几个阶段组成的过程链: ① 问题定义 ② 数据转换 ③ 数据探索 ④ 预测模型 ⑤ ...
  • 眼动数据分析基础_数据处理

    千次阅读 2020-12-21 19:10:25
    在记录眼动数据的时候,眼动仪按采样率来采集眼动原始数据(采样率30、60、120货300HZ)。每个数据点都将被识别为一个时间标签或者(x,y)的坐标形式并且被发送到与眼动以连接的加算机上运行的分析软件程序的数据库中...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 345,980
精华内容 138,392
关键字:

原始数据分析方法