0%

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
#####################################################################

写这个原因呢,最近又要对样本的HLA分子进行分型,然后看到某公司的微信公众号讲的HLA的分型软件,全文讲了那么多,要么巨难用,要么下载不到,反正不如我自己正在用的这两个。另外一方面,没必要太纠结非常高的精度,除非你用得到。4-digital resolution,我觉得已经够了。

HLA分子

先回顾下百度百科对HLA的介绍(https://baike.baidu.com/item/HLA/9504270?fr=aladdin):

HLA(human leukocyte antigen ,人类白细胞抗原)是人类的主要组织相容性复合体(MHC)的表达产物,该系统是所知人体最复杂的多态系统。

HLA是具有高度多态性的同种异体抗原,其化学本质为一类糖蛋白,由一条α重链(被糖基化的)和一条β轻链非共价结合而成。其肽链的氨基端向外(约占整个分子的3/4),羧基端穿入细胞质,中间疏水部分在胞膜中。HLA按其分布和功能分为Ⅰ类抗原和Ⅱ类抗原。

HLA-I类分子:内源性抗原的递呈分子, HLA-Ⅱ类分子:外源性抗原的递呈分子

我想介绍的是seq2hla和optitype这两个软件

软件的大致原理,都是和HLA数据库中的数据进行比对,然后分型。我拿过7个样本的RNA数据用这两个软件进行比较,如下表,虽然个别分型不一致(后两个resolution不一致),但大部分结果还是相当一致的。

SAMPLE A1 A2 B1 B2 C1 C2
1 seq2HLA A*02:01 A*24:18 B*35:102 B*35:102 C*16:02 C*04:01
OptiType A*02:01 A*24:02 B*35:14 B*35:14 C*16:02 C*04:01
2 seq2HLA A*29:01 A*02:07 B*46:01 B*15:01 C*04:01 C*01:02
OptiType A*29:01 A*02:07 B*46:01 B*15:01 C*04:01 C*01:02
3 seq2HLA A*02:07 A*11:01 B*52:01 B*15:01 C*12:02 C*03:03
OptiType A*02:07 A*11:01 B*52:01 B*15:01 C*12:02 C*03:03
4 seq2HLA A*02:07 A*11:01 B*46:01 B*40:01 C*01:02 C*07:02
OptiType A*02:07 A*11:01 B*46:01 B*40:01 C*01:02 C*07:02
5 seq2HLA A*33:03 A*11:01 B*58:01 B*15:02 C*03:02 C*08:01
OptiType A*33:03 A*11:01 B*58:01 B*15:02 C*03:02 C*08:01
6 seq2HLA A*02:07 A*29:01 B*46:01 B*07:02 C*01:02 C*15:05
OptiType A*02:07 A*29:01 B*46:01 B*07:05 C*01:02 C*15:05
7 seq2HLA A*02:03 A*11:01 B*56:01 B*56:01 C*01:02 C*01:02
OptiType A*02:03 A*11:01 B*56:01 B*56:01 C*04:01 C*01:02

seq2hla

首先介绍的是seq2hla,这个其实一个命令就可以进行HLA的分型(很实用吧),但只基于RNA Seq的数据,我见有人用WES的数据来做,但不确定靠不靠谱。

官网: https://bitbucket.org/sebastian_boegel/seq2hla/wiki/Home

利用conda安装

1
2
conda create -n seq2hla
conda install -n seq2hla -c bioconda seq2hla

运行

1
2
seq2HLA -1  readfile1  -2  readfile2  -r  runname  [-p int] [-3 int]
# 提供R1和R2,现在已经支持gzip格式,-p是线程,-3是trim read末端的碱基数(因为read末端的碱基质量不高)

seq2hla的结果

下面的结果基于RNA-seq的数据得到的,在介绍optitype的时候,我用的是WES-seq的数据,可以看到虽然软件和数据类型不同,但两个软件的结果是一致的,综合前面表格中两种软件基于同一样本RNA-seq的结果,和两种软件基于同一样本不同数据类型的结果来看,都能说明两个软件都比较靠谱。

1
2
3
4
5
6
7
8
#Locus  Allele 1        Confidence      Allele 2        Confidence
A A*24:02 0.001426754 A*11:01 0.003463472
B B*38:02' 0.0002728549 B*46:01 0.0002107604
C C*01:02 0.01282529 C*07:02 0.01516768
#Locus Allele 1 Confidence Allele 2 Confidence
DQA no NA no NA
DQB DQB1*05:02 NA DQB1*05:02 NA
DRB DRB1*11:52 0.297521 DRB1*04:05' 0.42202

optitype

另外一个软件是optitype,这个软件呢,稍微多了几步,但绝对在让人接受的范围内,而且这个软件支持DNA数据和RNA数据。

官网:https://github.com/FRED-2/OptiType

利用conda安装

这个时候真的要安利下conad,有了conda之后,各种软件包的安装和管理,变得非常轻松。拿optitype来说,我第一次安装的时候,自己下载,编译,修改部分文件,还要自己安装razers,而有了conda之后,只需一个命令

1
2
conda create -n optitype
conda install -n optitype -c bioconda optitype

optitype需要先利用razers3先过滤一下reads,把和HLA相关的序列提出来,然后在运行optitype进行分型,从而加快分析速度。

第一步,Razers3

1
razers3 -i 95 -m 1 -dr 0 --thread-count 10 -o fished_1.bam /path/to/hla_reference_dna.fasta sample_1.fastq

我是用conda装的,但在conda的env目录下没找到hla的fasta文件,不过不用着急,官网上有,可以从 https://github.com/FRED-2/OptiType/tree/master/data 下载对应的HLA的DNA或RNA的数据。

如果是paired的read,则需要分别跑一次,得到fished_1.bam和fished_2.bam文件。

另外官网提示Razers3会把所有的fastq读到内容,所以计算服务器的内存最好大点。

第二步,将bam文件转成fastq

1
2
samtools bam2fq fished_1.bam > fished_1.fastq
samtools bam2fq fished_2.bam > fished_2.fastq

第三步,运行optitype

1
OptiTypePipeline.py -i fished_1.fastq fished_2.fastq --dna -v -o /path/to/outdir -p sample.optitype.dna

–dna需要和前面Razers3用的HLA的序列类型一致,如果前面用的RNA的数据,这里也是–rna,当然用RNA还是DNA取决于测序数据。

optitype的结果
这个结果和seq2hla基于RNA-seq的数据得到的HLA亚型是一致的。

1
2
        A1      A2      B1      B2      C1      C2      Reads   Objective
0 A*11:01 A*24:02 B*38:02 B*46:01 C*01:02 C*07:02 200.0 191.00000000000003

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

安捷伦的全外芯片捕获的目标区域的bed文件是可以下载的,下载网址

https://earray.chem.agilent.com/suredesign

这个网站我没保存过,毕竟用的频率不高,但每次想用的时候还要搜一下这个网站(Σ( ° △ °|||)︴)

注册登录之后,找到Find Design -> SureSelect DNA -> Agilent Catalog Designs,进一步筛选之后,可以下载对应的文件。

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