十一 10

一个基因在同一套基因组上面竟然有两个定位!

查了好久的bug,终于搞清楚问题所在了!因为要对基因进行reads计数,所以要拿到基因在基因组上面的染色体起始终止坐标,结果发现了个十分诡异的现象,很多基因有多个坐标,比如下面这个PTPRS 在hg38这个基因组版本,居然有两个定位,因为我是写程序格式化得到的坐标,所以我check了我的程序,http://www.biotrainee.com/thread-472-1-1.html  感兴趣的同学可以点开看看我的代码!

tmp
Continue reading

23

我挣大钱了?

最近跟一些志同道合的小伙伴们一起搭建了 生信技能树 的论坛,所以在社交上面加大了投入力度,认识了很多在生物信息学领域各个学习程度的同学,发现有些人问我的问题让哭笑不得,大意就是:我的生信菜鸟团博客有近四百篇文章,阅读10万+了,又是知乎上面的大V,现在又在建论坛,感觉生意很红火的样子,是不是挣了很多钱啊!

Continue reading

03

生信技能树论坛诞生啦!!!

在许多小伙伴的共同协作下,我们的第一个论坛-生信技能树,诞生啦!

论坛地址:http://www.biotrainee.com/forum.php

虽然大家都说论坛已经是过气的互联网产品了,但我对互联网行业懂的很少,其实当初做博客的时候就有人跟我说过类似的话,但我还是坚持做了,我觉得做得还挺成功的,所以我仍然决定坚持把这个论坛做下去。

博客有很多缺点,传播速度很慢,不利于检索分类文章,个人知识面有限,也没办法跟follower及时交流。而我们的论坛,就可以克服那几个缺点。 Continue reading

14

讨论-用高通量测序方法研究sepsis

估计很多小伙伴都没有听过sepsis,现在翻译成中文是脓毒病,很多人会把它与那个缺乏维他命C的败血症混淆,其实完全不一样,因为sepsis致死率非常高!

sepsis [ˈsepsɪs] ['sepsɪs]
n. 脓毒病; 脓毒疾;
[例句]This may be of value in the treatment of meningitis and sepsis.
这可能会在治疗脑膜炎和败血症上有一定价值。

Continue reading

14

我也想开个公司(下)

自从我写了那篇关于创业的想法的文章后,传送门:  我也想开个公司(上)  , 很多熟悉的朋友,还有不少陌生的朋友都给我来信,跟我讨论我的想法,尤其是几个海外的朋友特别热情,我们深度的讨论了创建自由职业者联盟的可行性,公司如何活下去,盈利点是什么,什么样的价值观才能铸就百年企业等各种话题,可能我们在创业这个领域都还是蹒跚学步的状态,但是大家的热心帮助让我很感动,这也坚定了我继续做知识传播者的想法。我这里简单分享一下我们讨论的几个公司战略方向,因为我还有四五年的博士生涯要度过,暂时无法全心全意的发展事业,如果有人看到了也想朝这个方向发展,我可以免费做咨询服务,我非常乐意看到有人能实现我的想法。 Continue reading

06

如果你希望我回答你的问题

最近有很多朋友咨询我关于生物信息学数据处理的各种问题,有通过QQ直接对话聊天的,或者在QQ群里at我的,或者在知乎上面给我发短信息的,还有给我的163邮箱发信的。怎么说呢,首先还是感谢大家对我的信任,愿意花时间来跟我交流生物信息学数据处理的相关技术,然后我要简单说明一下为什么有些时候我没有答复你,虽然可能对你来讲,我是没有礼貌或者是太傲气了,但是我在这个领域浸淫了这么久,虽然你愿意跟我交流,但是你们的很多问题对我来说要么是都是太小儿科了,简单的google就能解决,要么是太空泛了,我无从答起,甚至我也给不出正确答案,更多的是有些人压根不用心的提问,纯粹是耽误你我的时间,所以我觉得很有必要写这篇博客简单说明一下,什么情况下我会回答你的问题。(如果你的问题非常吸引人,下面你就不用看了,我一定会抢着回答你的!) Continue reading

04

跟师妹聊Exome-seq、ChIP-seq、RNA-seq之间的差异

最近学习CHIP-seq的分析流程,略有点心得,也跟以前掌握的WES和RNA-seq做了一些比较,趁跑步的时候跟师妹讨论了一下,正好师妹写了一篇博客来分享这个讨论结果,我也借此机会转载过来,分享给大家,算是借花献佛吧!师妹的博文是用markdown写作,我觉得大家应该直接看她的文章,写得条理清楚:Exome-seq、ChIP-seq、RNA-seq之间的差异 Continue reading

23

生物医疗大数据高峰论坛参会笔记(全)

呀,这是去年(2015)蹭的一个论坛,不记得是第几届了,反正是生物谷举办的,他们搞论坛已经成为一个产业了,非常挣钱的!我那时候还很认真的做了笔记,现在回过头来看看,他们好像讲的都很有道理,虽然我直到现在也用不上,不过我丝毫不担心。我一直拼命的学习各种知识,就是因为有着坚定的信念,所学的一切终将会有一天对我的人生有所帮助。

Continue reading

06

突变频谱探究mutation siganures

这也是对TCGA数据的深度挖掘,从而提出的一个统计学概念。文章研究了30种癌症,发现21种不同的mutation signature。如果理解了,就会发现这个其实蛮简单的,他们并不重新测序,只是拿已经有了的TCGA数据进行分析,而且居然是发表在nature上面!

研究了4,938,362 mutations from 7,042 cancers样本,突变频谱的概念只是针对于somatic 的mutation。一般是对癌症病人的肿瘤组织和癌旁组织配对测序,过滤得到的somatic mutation,一般一个样本也就几百个somatic 的mutation。

paper链接是:http://www.nature.com/nature/journal/v500/n7463/full/nature12477.html

从2013年提出到现在,已经有30种mutation siganures,在cosmic数据库有详细记录,更新见:http://cancer.sanger.ac.uk/cosmic/signatures
它的概念就是:根据突变上下文分成96类,然后每类突变的频率不一样画一个条形图,可视化展现。
mutation signature

Each signature is displayed according to the 96 substitution classification defined by the substitution class and sequence context immediately 3′ and 5′ to the mutated base. The probability bars for the six types of substitutions are displayed in different colours.
仔细看paper,还是蛮好理解的,自己写一个脚本就可以做这个分析了,前提是下载各个癌症的somatic mutation文件,一般是maf格式的,很多途径下载。
In principle, all classes of mutation (such as substitutions, indels, rearrangements) and any accessory mutation characteristic, for example, the sequence context of the mutation or the transcriptional strand on which it occurs, can be incorporated into the set of features by which a mutational signature is defined. In the first instance, we extracted mutational signatures using base substitutions and additionally included information on the sequence context of each mutation. Because there are six classes of base substitution—C>A, C>G, C>T, T>A, T>C, T>G (all substitutions are referred to by the pyrimidine of the mutated Watson–Crick base pair)—and as we incorporated information on the bases immediately 5′ and 3′ to each mutated base, there are 96 possible mutations in this classification. This 96 substitution classification is particularly useful for distinguishing mutational signatures that cause the same substitutions but in different sequence contexts.

很多癌症都发现了不止一种mutation signature,甚至高达6种,说明癌症之间差异还是蛮大的!
In most cancer classes at least two mutational signatures were observed, with a maximum of six in cancers of the liver, uterus and stomach. Although these differences may, in part, be attributable to differences in the power to extract signatures, it seems likely that some cancers have a more complex repertoire of mutational processes than others.
Most individual cancer genomes exhibit more than one mutational signature and many different combinations of signatures were observed
但是,我最后也没能绝对的界限是什么,因为总不能用肉眼来看每个突变频谱不一样吧?
The set of signatures will be updated in the future. This will include incorporating additional mutation types (e.g., indels, structural rearrangements, and localized hypermutation such as kataegis) and cancer samples. With more cancer genome sequences and the additional statistical power this will bring, new signatures may be found, the profiles of current signatures may be further refined, signatures may split into component signatures and signatures may be found in cancer types in which they are currently not detected.
分类会持续不断更新,随着更多的cancer type和样本加入,新的signature会被发现,现有的signature也可能会被重新定义,或者被分割成多个小的signature
23

我的生信资源网盘分享

我重新整理了一下我当初生信入门时候收集的资料,希望对后学者有帮助,如果你以前储存了我的分享,请删掉以前的, 重新转这个完整版,以后不会更新了,入门的资料其实很多都是一样的!

公开链接:http://pan.baidu.com/s/1jIvwRD8 ,没有密码,直接根据需求自行下载

Continue reading

23

我的优酷视频网盘分享

因为视频教程架构出了点问题,所有不会更新了。现在公布所有下载链接,大家不需要去优酷观看了,太不清晰了,而且还有广告!

我的自频道-优酷视频 

视频百度云盘链接:http://pan.baidu.com/s/1jIvwRD8

下面是以前写的大纲,虽然我不做视频了,但是大家可以按照下面的大纲来自学,希望能对你有所帮助!!!

Continue reading

12

bioconductor中文社区招募站长

我已经构建好bioconductor中文社区的雏形,大家可以进去看看!

写在前面

突然发现我的bioconductor.cn这个域名都快要过期了!

哈哈,才想起一年前的计划到现在还没开始实施,实在不像我的风格,可能是水平到了一定程度吧,很多初级工作不像以前那样事无巨细的把关了。正好,借这个机会找几个朋友帮我一起完成这个bioconductor中文社区计划!

 

Continue reading

10

THis is a test for markdown

A wonderful tool to write blog !

Hi,all. I just started to use github and found this magic way to write a blog.
so I do a simpler test here.
Don't be confused, there is nothing new for bioinformatic here.

And it very easy to use Markdown in wordpress,just download a plugin.
Below is some useful links to help you familar with markdown as soon as possible.

  • http://bhttp://mahua.jser.me/log.csdn.net/kaitiren/article/details/38513715
  • https://github.com/guodongxiaren/README/blob/master/README.md
  • http://wowubuntu.com/markdown/
  • http://mahua.jser.me/

Also I want to recommend some useful web-editor for markdown.

  • http://mahua.jser.me/
  • https://www.madoko.net/
  • https://www.zybuluo.com/mdeditor

It's easy to use Markdown, but to get the most out of it, you still need to understand it and keep practising

22

没必要学shell进阶语法

因为大部分生物信息学软件都是linux版本的,所以生物信息学数据分析工作者必备技能就是linux,但是大部分人只是拿他当个中转站,我以前也是,直到接触了大批量的任务,自动化流程,才明白这里面的水太深了,不过无所谓,凭我个人的观点,其实shell的进阶语法真的不必要!

当然,只是我一家之言!
我实在是不想去背诵大括号,小括号,中括号以及双重括号到底区别是什么!
http://www.bio-info-trainee.com/?p=1018  [],[[]],(),(()),{},{{}},以及在前面加上$的区别,以及它们互相杂交组合的区别!!!

我也不想去搞明白操作符两边是否加空格的区别是什么了。

if((i%5==0)) 来判断变量是否被一个数整除
i=$((i+1))来表示变量自增。
这些东西真的很诡异!
如果你有qsub,condor等任务提交系统,那么你只需要熟悉他们就可以了,但大部分散兵游勇的生物信息学家并没有集群,所以压根不会接触任务提交系统,就需要些自动化脚本了!
受限制与机器的cpu以及内存数,需要判断提交了多少任务,等待多久再执行,所以会把一个简单的自动化脚本写的很复杂!
比如下面这个脚本:cat >download_hg38_from_UCSC.sh

for i in $(seq 1 22) X Y M;
do echo $i;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg38.fa;
done
rm -fr chr*.fa
可以下载hg38基因组的fasta文件,但是是分染色体一个个下载的!
再比如下面这个,批量做GSEA分析的脚本:
while read id
do
echo $id
gene=`echo $id |awk '{print $1}'`
probe=`echo $id |awk '{print $2}'`
echo $i
do_GSEA $probe $gene; ##这里是我自己定义的一个function,就不贴出来了
if((i%5==0))
then
sleep 10  ##重点就在这里,每次提交的任务有限制,所以需要休息,不然机器的cpu负载太高!
fi
i=$((i+1))
done <$1

如果,还有其它功能需要实现,我们可以把脚本写的更负载,纯粹的用shell,需要搜索更多的shell技巧。

但是事实上并没有这个必要,我们现在有了更方便的脚本语言,比如我所擅长的perl
我写一个nohup提交任务的脚本!
## perl nohup.pl   deep_count.sh  0
## perl nohup.pl   deep_count.sh  1
## perl nohup.pl   deep_count.sh  2

[perl]
## perl nohup.pl   deep_count.sh  0
## perl nohup.pl   deep_count.sh  1
## perl nohup.pl   deep_count.sh  2
$i=1;
open FH,$ARGV[0];
while(<FH>){
    chomp;
    next unless $.%3==$ARGV[1];
    $cmd="nohup  $_  &";
    print "$cmd\n";
    system($cmd);
    sleep(10800) if $i%5==4;
    $i++;
    #exit;
}
[/perl]

我尝试过用shell,写了很久,总是报错,但是用perl,一分钟我就写完了,所以,最好是用自己熟悉的一种语法最好!
08

别人写的代码运行真快!!!

最近需要做十几万个差异基因分析,每个分析都是对大约5万个探针,200个样本的数据量进行批量T检验计算P值
然后发现自己无论怎么用R来写,每个分析都要耗时半分钟左右,因为我必须循环所有的探针,即使不用for,而用R推荐的apply系列函数,也快不到哪里去,但是我搜索时候发现有一个package里面自带了矩阵T检验,直接对5万个探针进行T检验,而不需要循环处理它们
看下面代码!
dat=matrix(rnorm(10000000),nrow = 50000)
dim(dat) #50000   200
system.time(
  apply(dat,1,function(x){
    t.test(x[1:100],x[101:200])$p.value
  })
)
#用户  系统  流逝
#29.29  0.04 30.64
library(pi0)
system.time(matrix.t.test(dat,1,100,100))
#用户 系统 流逝
#0.48 0.03 0.53
差距真的是非常的明显呀!!!
然后,我解析了它的代码,发现里面调用了C写的代码,我想这就是问题所在咯,可是他们到底怎么写,才能把速度搞这么快???
  tmp = .C("tstatistic", dat = x, n1 = n1, n2 = n2, ntests = ntests, 
        MARGIN = MARGIN, pool = pool, tstat = rep(0, ntests), 
        df = rep(0, ntests), PACKAGE = "pi0")

源码在这个package的github里面可以找到,有兴趣的童鞋可以研究一下

 
 

 

十二 15

用超几何分布检验做富集分析

我们可以直接使用R的bioconductor里面的一个包,GOstats里面的函数来做超几何分布检验,看看每条pathway是否会富集

我们直接读取用limma包做好的差异分析结果

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

DEG=read.table("GSE63067.diffexp.NASH-normal.txt",stringsAsFactors = F)

View(DEG)

image001

我们挑选logFC的绝对值大于0.5,并且P-value小雨0.05的基因作为差异基因,并且转换成entrezID

probeset=rownames(DEG[abs(DEG[,1])>0.5 & DEG[,4]<0.05,])

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

length(unique(EGID))

#[1] 775

diff_gene_list <- unique(EGID)

这样我们的到来775个差异基因的一个list

首先我们直接使用R的bioconductor里面的一个包,GOstats里面的函数来做超几何分布检验,看看每条pathway是否会富集

library(GOstats)

library(org.Hs.eg.db)

#then do kegg pathway enrichment !

hyperG.params = new("KEGGHyperGParams", geneIds=diff_gene_list, universeGeneIds=NULL, annotation="org.Hs.eg.db",

categoryName="KEGG", pvalueCutoff=1, testDirection = "over")

KEGG.hyperG.results = hyperGTest(hyperG.params);

htmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))

结果如下:

image002

但是这样我们就忽略了其中的原理,我们不知道这些数据是如何算出来的,只是由别人写好的包得到了结果罢了。

事实上,这个包的这个hyperGTest函数无法就是包装了一个超几何分布检验而已。

如果我们了解了其中的统计学原理,我们完全可以写成一个自建的函数来实现同样的功能。

超几何分布很简单,球分成黑白两色,数量已知,那么你随机抽有限个球,应该抽多少白球的问题!
公式就是 exp_count=n*M/N
然后你实际上抽了多少白球,就可以计算一个概率值!
换算成通路的富集概念就是,总共有多少基因,你的通路有多少基因,你的通路被抽中了多少基因(在差异基因里面属于你的通路的基因),这样的数据就足够算出上面表格里面所有的数据啦!
tmp=toTable(org.Hs.egPATH)
GeneID2Path=tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)
Path2GeneID=tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
#phyper(k-1,M, N-M, n, lower.tail=F)
#n*M/N
diff_gene_has_path=intersect(diff_gene_list,names(GeneID2Path))
n=length(diff_gene_has_path) #321 # 这里算出你总共抽取了多少个球
N=length(GeneID2Path) #5870  ##这里算出你总共有多少个球(这里是错的,有多少个球取决于背景基因!一般是两万个)
options(digits = 4)
for (i in names(Path2GeneID)){
 M=length(Path2GeneID[[i]])  ##这个算出你的所有的球里面,白球有多少个
 exp_count=n*M/N  ###这里算出你抽取的球里面应该多多少个是白色
 k=0         ##这个k是你实际上抽取了多少个白球
 for (j in diff_gene_has_path){
 if (i %in% GeneID2Path[[j]]) k=k+1
 }
 OddsRatio=k/exp_count
 p=phyper(k-1,M, N-M, n, lower.tail=F)  ##根据你实际上抽取的白球个数,就能算出富集概率啦!
 print(paste(i,p,OddsRatio,exp_count,k,M,sep="    "))
}
随便检查一下,就知道结果是一模一样的!

 

十二 13

可怕的中国教育(转)

无意中看到这篇文章,感觉挺有趣的!
我不知道国外的教育是咋样的,但是文中描述的中国教育的试题都是千真万确的,我也经历过。
我没看过有钱人家的小孩,不知道那些花了几百万买个几平米的学区房的那些人的小孩是不是也要经受这样的教育摧残?
我还有个问题,我是如何出淤泥而不染的呢?
随便瞎扯,大家自己看这篇文章吧!
可怕的中国教育,聪明伶俐进去呆若木鸡出来
       作者:葛小琴
  来源:杂文月刊(文摘版)
  侄子在读高二,考了一道历史题:成吉思汗的继承人窝阔台,公元哪一年死?最远打到哪里?答不出来,我帮他查找资料,所以到现在我都记得,是打到现在的匈牙利附近。
  在一次偶然的机会,我发现美国世界史这道题目不是这样考的。它的题目是这样的:成吉思汗的继承人窝阔台,当初如果没有死,欧洲会发生什么变化?试从经济、政治、社会三方面分析。
   有个学生是这样回答的:这位蒙古领导人如果当初没有死,那么可怕的黑死病,就不会被带到欧洲去,后来才知道那个东西是老鼠身上的跳蚤引起的鼠疫。但是六 百多年前,黑死病在欧洲猖獗的时候,谁晓得这个叫做鼠疫?如果没有黑死病,神父跟修女就不会死亡。神父跟修女如果没有死亡,就不会怀疑上帝的存在。如果没 有怀疑上帝的存在,就不会有意大利弗罗伦斯的文艺复兴?
  如果没有文艺复兴,西班牙、南欧就不会强大,西班牙无敌舰队就不可能建立。如果西班牙、意大利不够强大,盎格鲁—撒克逊会提早200年强大,日耳曼会控制中欧,奥匈帝国就不可能存在。
  教师一看“棒,分析得好。”但他们没有分数,只有等级A。其实这种题目老师是没有标准答案的,可是大家都要思考。
  不久前,我去了趟日本,日本总是和我们在历史问题上产生纠葛,所以我在日本很注意高中生的教科书。
   他们的教师给高中生布置了这样一道题:日本跟中国100年打一次仗,19世纪打了日清战争(即甲午战争),20世纪打了一场日中战争(即抗日战 争),21世纪如果日本跟中国开火,你认为大概是什么时候?可能的远因和近因在哪里?如果日本赢了,是赢在什么地方?输了是输在什么条件上?分析之。
   其中有个高中生是这样分析的:我们跟中国很可能在台湾回到中国以后,有一场激战。台湾如果回到中国,中国会把基隆与高雄封锁,台湾海峡就会变成中国的内 海,我们的油轮就统统走右边,走基隆和高雄的右边。这样,会增加日本的运油成本。我们的石油从波斯湾出来跨过印度洋,穿过马六甲海峡,上中国南海,跨台湾 海峡进东海到日本海,这是石油生命线,中国政府如果把台湾海峡封锁起来,我们的货轮一定要从那里经过,我们的主力舰和驱逐舰就会出动,中国海军一看到日本 出兵,马上就会上场,就开打。
  按照判断,公元2015年至2020年之间,这场战争可能爆发。所以,我们现在就要做对华抗战的准备。
  我看其他学生的判断,也都是中国跟日本的摩擦会从东海从台湾海峡开始,时间判断是2015年至2020年之间。
  这种题目和答案都太可怕了。
   撇开政治因素来看这道题,我们的历史教育就很有问题。翻开我们的教科书,题目是这样出的:甲午战争是哪一年爆发的?签订的叫什么条约?割让多少土地?赔 偿多少银两?每个学生都努力做答案。结果我们一天到晚研究什么时候割让辽东半岛,什么时候丢了台湾、澎湖、赔偿二万银两,1894年爆发甲午战争、 1895年签订马关条约,背得滚瓜烂熟,都是一大堆枯燥无味的数字。
  那又怎么样,反正都赔了嘛,银两都给了嘛,最主要的是将来可能会怎样。
  人家是在培养能力,而我们是在灌输知识,这是值得深思的部分!
  看外面的教育,再看我们的教育?
  老妈去参加我侄子的家长会,带回了两套侄子的考试试卷,我很好奇,拿过来看了现在小学生的试卷后,我震惊了!这是什么狗屁教育?这样的教育有希望吗?下面给大家详细说说我看到了什么。
  侄子在本市某著名小学读书,有这么几道题。
  一个春天的夜晚,一个久别家乡的人,望着皎洁的月光不禁思念起了故乡,于是吟起了一首诗:( ),( )。
   我看到侄子答的是:举头望明月,低头思故乡。但后面是一把大大的×,我就奇怪了,我也是想到的这两句。好奇地问侄子,这个不对??那答案是什么?侄子说 标准答案是:春风又绿江南岸,明月何时照我还。哎,这就奇怪了,因为是个春天的夜晚,就要是这句有春风的???要这个思念故乡的人不是江南的,是不可能说 出春风又绿江南岸这句话的!!!举头望明月,低头思故乡应该更准确。再扯远点,思念故乡,一千个人可以吟一千句不一样的诗,这个也可以有标准答案的么?
  接下来是默写,题目是:我们学过《桂林山水》一文,请将下面句子默写下来。然后就是整段的要默写,这有什么用?死记硬背别人的文字有什么用?
   还有个题目,《匆匆》这篇课文,是现代着名作家朱自清先生写的,同学们都很喜欢这篇散文,你能把自己最喜欢,印象最深刻的一句写下来吗?我侄子写的是: 我的日子滴在时间的流里,没有声音,也没有影子。后面一把好大的×。标准答案竟然是:但是,聪明的你告诉我,我们的日子为什么一去不复返呢?这就更奇怪 了,一篇文章,你可以喜欢这句,我可以喜欢那句,难道最喜欢的一句话也要统一么?为什么“我的日子滴在时间的流里,没有声音,也没有影子。”这句不能喜 欢?就一定要喜欢“但是,聪明的你告诉我,我们的日子为什么一去不复返呢?”这句?我觉得这个题目应该是“你能把老师最喜欢,印象最深刻的一句写下来 吗?”才对!
  再看别的试卷,更莫名其妙了,比如请说出阿拉伯数字的来历,是哪个国家创造的?侄子不知道,问我,我也不知道。我只好去搜一下,才知道是古印度人发明的。莫非我吃块猪肉,还一定得知道它是哪个养猪场养出来的?
   最后有个题目让我彻底崩溃了:请用一句话说明π的含义。侄子回答π的含义是圆周率。竟然打的是×,这就奇怪了,正好我老婆大学说读的是理科,我马上问 她,π是什么意思,她说圆周率啊。两个人狂汗,问了侄子半天,标准答案大概是,π是一个在数学及物理学领域普遍存在的数学常数……
  如果你也觉得这种教育很无耻,就请转发吧,让更多的人来参与呼吁改变,为了孩子为了国家的未来 …… 别让孩子聪明伶俐地进去呆若木鸡地出来!
十二 11

关于limma包差异分析结果的logFC解释

首先,我们要明白,limma接受的输入参数就是一个表达矩阵,而且是log后的表达矩阵(以2为底)。

那么最后计算得到的logFC这一列的值,其实就是输入的表达矩阵中case一组的平均表达量减去control一组的平均表达量的值,那么就会有正负之分,代表了case相当于control组来说,该基因是上调还是下调。

我之前总是有疑问,明明是case一组的平均表达量和control一组的平均表达量差值呀,跟log foldchange没有什么关系呀。

后来,我终于想通了,因为我们输入的是log后的表达矩阵,那么case一组的平均表达量和control一组的平均表达量都是log了的,那么它们的差值其实就是log的foldchange

首先,我们要理解foldchange的意义,如果case是平均表达量是8,control是2,那么foldchange就是4,logFC就是2咯

那么在limma包里面,输入的时候case的平均表达量被log后是3,control是1,那么差值是2,就是说logFC就是2。

这不是巧合,只是一个很简单的数学公式log(x/y)=log(x)-log(y)