0%

datapasta

https://github.com/MilesMcBain/datapasta/

还在手工的把excel的数据写成导到R里吗。不管横着还是竖着复制数据,datapasta可以自动、快速的把复制数据转成tibbles, data.frames, 或者 vectors格式。

更详细的参考https://cran.r-project.org/web/packages/datapasta/vignettes/how-to-datapasta.html

https://raw.githubusercontent.com/lorenzwalthert/some_raw_data/master/styler_0.1.gif

styler

https://github.com/r-lib/styler

还嫌自己写出的代码不够美吗,styler 可以把代码格式化成tidyverse规则的风格。

The goal of styler is to provide non-invasive pretty-printing of R source code while adhering to the tidyverse formatting rules. styler can be customized to format code according to other style guides too.

https://raw.githubusercontent.com/lorenzwalthert/some_raw_data/master/styler_0.1.gif

真香

links

https://github.com/MilesMcBain/datapasta/

https://github.com/r-lib/styler

https://github.com/daattali/addinslist

https://www.zhihu.com/question/398418315

新年快乐,21年的第一篇文章。

以前写过映射ENSEMBL ID 和 NCBI ID, http://www.zxzyl.com/archives/736

日常分析中,我们也会经常遇到其他的ID mapping的工作,这种工作不是基因ID转基因ID,而是转录本的ID转基因ID。

如果用的是refGene的注释,最简单了,直接用下面的命令即可

1
mysql --user=genome -N --host=genome-mysql.cse.ucsc.edu -A -D hg38 -e "select name,name2 from refGene"

不过我也经常通过解析gtf文件获得,因为gtf有转本的ID,也有基因的symbol或者ID,只要有gtf文件就可以提取。本着不造轮子的精神,我利用的是现成的R包

1
2
3
library(plyranges)
gr <- read_gff("/path/to/gtf/or/gff") %>% select(transcript_id, gene_id, gene_name)
gr <- unique(data.frame(gr))

我也在自己的包里面写了一个函数得到ensembl,refseq,hgnc,gene symbol的对应关系,biomaRt比较慢,可以把结果保存成文件

1
2
3
4
5
devtools::install_github("ProfessionalFarmer/loonR")
# 需要安装biomaRt和dplyr
mapping.table <- loonR::get_full_mapping_table()
# 保存mapping.table

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

Fixation index (FST)

整理来源 http://www.uwyo.edu/dbmcd/popecol/maylects/fst.html,把这个课程的步骤用表格来表示

Subpopulation 1 Subpopulation 2 Subpopulation 3 Total
Genotype AA 125 50 100
Aa 250 30 500
aa 125 20 400
Number of individual 500 100 1000 1600
Number of alleles 1000 200 2000 3200
Step 1. Calculate the gene (allele) frequencies
Observed allele frequency A (p) (125*2+250)/1000=0.5 (50*2+30)/200=0.65 (2*100+500)/2000=0.35
a (q) 0.5 0.35 0.65
Step 2. Calculate the expected genotypic counts under Hardy-Weinberg Equilibrium, and then calculate the excess or deficiency of homozygotes in each subpopulation. Summary of homozygote deficiency or excess relative to HWE: Pop. 1. Observed = Expected: perfect fit Pop. 2. Excess of 15.5 homozygotes: some inbreeding Pop. 3. Deficiency of 45 homozygotes: outbred or experiencing a Wahlund effect (isolate breaking).
Expected allele frequency AA 500*0.5^2 = 125 (= observed) 100*0.65^2 = 42.25 (observed has excess of 7.75) 1,000*0.35^2 = 122.5 (observed has deficiency of 22.5)
Aa 50020.5*0.5 = 250 (= observed) 10020.65*0.35 = 45.5 (observed has deficit of 15.5) 1,00020.65*0.35 = 455 (observed has excess of 45)
aa 500*0.5^2 = 125 (= observed) 100*0.35^2 = 12.25 (observed has excess of 7.75) 1,000*0.35^2 = 422.5 (observed has deficiency of 22.5)
Step 3. Calculate the local observed heterozygosities of each subpopulation (we will call them Hobs s, where the s subscript refers to the sth of n populations – 3 in this example).
Local observed heterozygosities 250/500 = 0.5 (Hobs 1) 30/100 = 0.3 (Hobs 2) 500/1000 = 0.5(Hobs 3)
Step 4. Calculate the local expected heterozygosity, or gene diversity, of each subpopulation Hexp = 2pq
Local expected heterozygosity 20.50.5=0.5 (Hexp 1) 20.653.5=0.455 (Hexp 2) 20.350.65=0.455 (Hexp 2)
Step 5. Calculate the local inbreeding coefficient of each subpopulation F = (Hexps -Hobs)/Hexp [positive F means fewer heterozygotes than expected indicates inbreeding] [negative F means more heterozygotes than expected means excess outbreeding]
F1=(0.5—0.5)/0.5=0 F2=(0.455—0.3)/0.455=0.341 F3=(0.455—0.5)/0.455=-0.099
Step 6. and 7. Calculate p-bar (p-bar, the frequency of allele A) over the total population. Calculate q-bar (q-bar, the frequency of allele a) over the total population. Check: p-bar + q-bar = 1.0
the frequency of allele over the total population p-bar (0.51000+0.65200+0.35*2000)/3200=0.4156
q-bar (0.51000+0.35200+0.65*2000)/3200=0.5844
Step 8. Calculate the global heterozygosity indices (over Individuals, Subpopulations and Total population) HI based on observed heterozygosities in individuals in subpopulations HS based on expected heterozygosities in subpopulations HT based on expected heterozygosities for overall total population
HI (observed) (0.5500+0.3100+0.5*1000)/1600=0.4875
HS (expected) (0.5500+0.455100+0.455*1000)/1600=0.4691
HT (in overall total population) 2*p-bar *q-bar = 2 * 0.4156 * 0.5844 = 0.4858
Step 9. Calculate the global F-statistics Compare and contrast the global FISbelow with the “local inbreeding coefficient” Fs of Step 5. Here we are using a weighted average of the individual heterozygosities over all the subpopulations. Both FIS and Fs are, however, based on the observed heterozygosities, whereas FST and FIT are based on expected heterozygosities.
FIS (Hs-Hi)/Hs=(0.4691-0.4875)/0.4691=-0.0393
FST (Ht-Hs)/Ht=(0.4858-0.4691)/0.4858=-0.0344
FIT (Ht-Hi)/Ht=(0.4858-0.4875)/0.4858=-0.0036
Step 10 conclusions

Finally, draw some conclusions about the genetic structure of the population and its subpopulations. 1) One of the possible HWE conclusions we could make: Pop. 1 is consistent with HWE (results of Step

  1. Two of the possible “local inbreeding” conclusions we could make from Step 5:

Pop. 2 is inbred (results of Step 5), and

Pop. 3 may have disassortative mating or be experiencing a Wahlund effect (more heterozygotes than expected).

  1. Conclusion concerning overall degree of genetic differentiation (FST)

Subdivision of populations, possibly due to genetic drift, accounts for approx. 3.4% of the total genetic variation (result of Eqn FST.8 FST calculation in Step 9),

  1. No excess or deficiency of heterozygotes over the total population (FIT is nearly zero).

I always use gtf file and retrieve gene information. There isn’t a highly flexible tool to solve my demand. I modified the code from “https://github.com/Jverma/GFF-Parser”, thanks Jverma. This tool will be easier to use.

Usage

Basically, there are three parameters.

id: either transcript id or gene id.

attType: attribute defined in gtf file. E.g. feature (column 3), gene_name, exon_number, transcript_id in column 9

attValue: the attribute value you want to search for.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
>>> import sys
>>> from gtfParser import gtfParser

>>> gtf = gtfParser("example.gtf")

>>> # Get all exons in CDK7
>>> gtf.getRecordsByID("CDK7", "feature", "exon")

>>> # Get all features of transcript_id defined as "NM_001324069" in gene "CDK7"
>>> gtf.getRecordsByID("CDK7", "transcript_id", "NM_001324069")

>>> # Get start codon where feature was defined as "start_codon" in transcript "NM_001324069"
>>> gtf.getRecordsByID("NM_001324069", "feature", "start_codon")

>>> # Get a exon where its id is "NM_001324078.1" in "NM_001324078" transcript
>>> gtf.getRecordsByID("NM_001324078", "exon_id", "NM_001324078.1")
Read more »

我是在这篇文章(Integrative pathway enrichment analysis of multivariate omics data
)中遇到的合并多个p-value的操作。这篇文章是今年发表在NC上。所有的组学或者大规模的数据分析,都需要探索数据背后相关的生物学功能,所以通路富集分析非常普遍。通常的做法是基于单一组学、单一数据集的数据进行分析,随着生物学数据的爆发,大规模多组学数据变得普遍,这篇文章介绍了基于整合的多组学或多数据集的数据进行通路分析的工具ActivePathways。

方法

ActivePathways的方法,如下图:

(a) 需要的输入文件

(1) 基于多组学数据集的基因P-value,传统的富集分析是单组学,只有一列,现在是多组学,对应多列P-value
(2) 基因集,这个和其他的通路富集分析一样,用来表示生物学过程和通路

(b)

(1) 用Brown method合并基因的P-values,并且排序,用一个宽松的阈值来过滤检阳性的基因。
(2) 对每个通路,用排序的基因(从第一个开始从少到多作为sub-list)进行超几何检验,并找到最优的sub-list长度。
(3) 基因单一组学的数据进行富集分析,找到支持每个通路的证据。

(c) ActivePathways 提供整合之后的富集分析结果,相关的Brown P-value,支持通路的证据。还可以在Cytoscape中画Enrichmentmap的图,来分析更广泛的生物学主题。点为通路,边表示有共有基因。

例子

自然发到了NC这种水平的期刊上,出了豪华的团队外,ActivePathways一定有很厉害的功能才行。这篇文章提了好几个case study,我挑了一个,稍微讲一下,如下图。

昨天分析了乳腺癌中与预后相关的通路和过程,其中整合了三种数据集,TC(Tumor cell mRNA),TAC(Tumor adjacent cell mRNA)和CNA(拷贝数变异),这里也挺有意思,把TAC也纳入到P-values的矩阵中,并不是三个组学。

(a) Enrichment map,其中蓝色的表示仅由整合数据发现的通路,比如凋亡等,显示出其强大的功能。
(b) 乳腺癌Basel和HER2亚型中,与预后相关的免疫基因的Hazard Ratio(HR),显示出两种亚型间不同的免疫pattern。
(c) HER2亚型中,与预后相关的基因(凋亡过程负调控)在每个数据集上的P-value和整合后的Brown P-value。其中DUSP1在三个数据集中都非常显著,如果作者聚焦到了这个基因。
(d) 这个基因在肿瘤细胞mRNA,癌旁细胞mRNA和拷贝数变异三种数据集中,分成high和low两组,做生存分析,可以发现DUSP1低表达,预后显著好于高表达,以此证明ActivePathway的强大。

ActivePaths的下载地址:https://github.com/reimandlab/ActivePathways。为一个R包,整理好P-values的数据框之后,一步命令即可分析,此外结果还可以在Cytoscape上用Enrichment map展示。

1
2
# scores是P-values的数据框,GMT是基因集
ActivePathways(scores, fname_GMT)

参考

https://www.nature.com/articles/s41467-019-13983-9

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

今天在看文章的时候,发现原来p-value也可以合并。比如一个基因在不同组学数据的检验中对应了多个p-value,可以合并成一个。

常用的是Fisher’s method,

![](https://raw.githubusercontent.com/ProfessionalFarmer/f4w/master/2020/2020-09-11-Fisher method.svg)

-2[ln(P1) + ln(P2) + … + ln(Pi)]符合X2分布(自由度为2k,k为p-value的个数)。

还有Brown’s methods和 Kost’s methods,具体的介绍如下图。

![](https://raw.githubusercontent.com/ProfessionalFarmer/f4w/master/2020/2020-09-11-Combining dependent P-values.png)

参考:

https://academic.oup.com/bioinformatics/article/32/17/i430/2450768

https://www.nature.com/articles/s41467-019-13983-9

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

M值和B值的计算公式

https://link.springer.com/article/10.1186/s41241-017-0041-9

The relationship curve between M-value and Beta-value

M值和B值的对应关系

The histograms of Beta-value (left) and M-value (right) (27578 interrogated CpG sites in total)

M值和B值的分布

minfi包有getM和getBeta来分别计算M-values和Beta-values,包的作者认为,

  • M-values具有更好的统计特性,更适合用于进行下游的统计分析(差异分析等)
  • Beta-values更加容易解释,更能说明生物学上的意义
  • CHAMP包在load的时候,可以指定计算Beta-value还是M-value

一般来说,具体的β值的意义是:

  • 任何等于或大于0.6的β值都被认为是完全甲基化的
  • 任何等于或小于0.2的β值被认为是完全未甲基化的
  • β值在0.2和0.6之间被认为是部分甲基化的。

参考:

https://zhuanlan.zhihu.com/p/108364645

https://link.springer.com/article/10.1186/s41241-017-0041-9

我们的存储服务器有两组RAID,容量均大于150T,我在mount的时候,提示我

1
2
3
4
5
NTFS signature is missing.
Failed to mount '/dev/sdc': Invalid argument
The device '/dev/sdc' doesn't seem to have a valid NTFS.
Maybe the wrong device is used? Or the whole disk instead of a
partition (e.g. /dev/sda, not /dev/sda1)? Or the other way around?

是因为没有分区导致的,分区之后就可以了。分区的命令

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# 使用parted命令进行分区,等同parted; select /dev/sdc
parted /dev/sdc

# 创建分区表
mklabel gpt

# 使用print命令查看当前分区情况
print

# 留1M的空余空间,目的是为了让数据块整齐,提高磁盘的运行效率, -1表示分区的结尾 意思是划分整个硬盘空间为主分区
mkpart primary 1 -1

p # print的简写

# 使用q命令退出,
quit

# 退出之后会提示
会提示Information: You may need to update /etc/fstab.


# 格式化分区,为分区写入文件系统,格式为ext4
mkfs –t ext4 /dev/sdc1 # 格式化分区

# 使用blkid命令,找到 UUID,然后编辑 /etc/fstab,实现自动挂载
vim /etc/fstab

UUID=****** directory ext4 defaults 0 0

参考:

https://www.cnblogs.com/kreo/p/9462641.html

https://www.cnblogs.com/saszhuqing/p/9964262.html

确定物理网口对应的名称

在一台ubuntu的机器上,有四个物理网口,我想知道每个网口对应的MAC地址。使用ip a可以看到网口的MAC地址和名称,比如列出了ens1f0, ens1f1, ens4f0, ens4f1。
原来的网卡interface都是eth开头,后来改成了enp, ens等。

Names incorporating Firmware/BIOS provided index numbers for on-board devices (example: eno1)
Names incorporating Firmware/BIOS provided PCI Express hotplug slot index numbers (example: ens1)
Names incorporating physical/geographical location of the connector of the hardware (example: enp2s0)
Names incorporating the interfaces’s MAC address (example: enx78e7d1ea46da)
Classic, unpredictable kernel-native ethX naming (example: eth0)

那么如何确定机器上的ens1f0对应的哪个物理网口呢,可以用ethtool来实现,ethtool是用于查询及设置网卡参数的命令。用ethtool -p enos1f1,看哪个网口在闪灯,就能确定这个物理网口对应的名称。记得不要插网线。

1
ethtool -p|--identify DEVNAME   Show visible port identification (e.g. blinking)

如果没有一个网口亮灯,很可能是因为网口不支持,则可以尝试ethtool -t enosf1f1,大概在4秒之后,网口的灯会亮,这个时候就可以确定enos1f1对应的具体的物理网口了。

1
ethtool -t|--test DEVNAME       Execute adapter self test

很简单的一个命令,知道了就很简单,不知道就很难想到。

配置静态IP

1
2
3
4
5
6
7
8
9
10
/etc/netplan/00-installer-config.yaml

network:
ethernets:
enp0s3: # 网卡名
addresses: [192.168.1.3/24] # ip地址和子网掩码,24对应255.255.255.0
gateway4: 192.168.1.1 # 网关
nameservers:
addresses: [4.2.2.2, 8.8.8.8] # DNS
version: 2

配置好了之后,生效

1
sudo netplan apply

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

在SNV分析中,我们在算signature和样本mutation spectrum之间的相似性时,会用到cosine similarity。cosine similarity (distance)的公式,其实就是两个向量的夹角的cosine值,计算公式如下

它与欧式距离的差别如下图,cosθ就是similarity,而d则是欧氏距离Euclidean distance。

有些时候,距离也算作一种相似性,因为距离越远,说明两个样本越不相似。Euclidean distance和cosine similarity要根据情况来选择,最重要的是,是否要考虑weight or magnitude,参考下图。

在文本挖掘分析的时候,计算两个文本的相似性,我们可以统计每个词出现的次数,然后计算相似性(距离),因为文章有短有长,如果考虑单词出现的次数,那么字数多的文章一定与字数少的文章不一样(欧氏距离),所以如果我们不考虑这个量(magnitude)的时候,用cosine计算更加合适,结果也与欧氏距离不一样。

参考:

https://cmry.github.io/notes/euclidean-v-cosine#:~:text=While%20cosine%20looks%20at%20the,to%20actually%20measure%20the%20distance.

https://www.machinelearningplus.com/nlp/cosine-similarity/

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################