shannon指数的python实现
这里以稀释曲线的方式求香农指数。可以查看结果是否到达平台期来判断香农指数是否达到最大。
先来看看shannon指数的公式。
\(H = -\sum_{i=i}^{s}{p_i}{log_2}{p_i}\)
这里s是物种OTU，pi是该OTU的丰度。一般使用log2，也有使用ln的。
我的写法是这样的
# 首先将所有物种OTU按照出现次数（即reads数）加入到一个list中
# 比方说就是OTU1如果有10000条reads支持，在list中就有10000个OTU1元素
# 以下以OTU_list表示这个list
import random
from collections import Counter
import math
selectNum = 0
while selec...
从OTU矩阵到Observed_species稀释曲线
这里是Observed species 稀释曲线的获得方法。它是利用已测得16S rDNA序列中已知的各种OTU的相对比例，来计算抽取n个（n小于测得reads序列总数）reads时出现OTU数量的期望值，然后根据一组n值（一般为一组小于总序列数的等差数列）与其相对应的OTU数量的期望值做出曲线来。
从这个原理看，可以知道的是我们只要定一个随机梯度，对OTU进行随机选取就可以做出这个Observed species的稀释曲线矩阵了。
一般的，OTU表格长这样
首先使用python将这个表格处理成用来作图的矩阵
import random
OTU = open("OTU.txt", "r")
ob = open("observed_species.txt", "w")
ob.w...
weka（四舍五入就算是会机器学习了吧）
WEKA
Weka是一个University of Waikato用java编写的开源机器学习软件，有着自己的GUI界面，内置了很多算法工具，适合我这种半吊子手残。IBM有一篇比较详细的教程了。
在官网下载好weka，接下来使用经典的鸢尾花卉数据集进行测试。这个数据集有150条数据，几个数据分别是sepal_length，sepal_width，petal_length，petal_width以及class。接下来，随机选130条作为训练集，剩下的20条作为测试集，测试一下对class的预测。
数据处理
然后需要将数据处理成arff格式，像下面这样，这里最好保证class在第一列或者最后一列，方便后面使用。
@RELATION iris
@ATTRIBUTE sepal_le...
多序列比对
以前上课的时候，老师教我们用的是EMBL的在线的Clustal Omega，现在试试自己搭一套。
多序列比对的软件有Clustal Omega、Kalign、MAFFT、MUSCLE、T-coffee、PRANK等等。
软件
这里选择使用Clustal Omega，然后再用mview进行可视化。
clustalo
Clustal Omega在页面上选择自己的版本下载就可以了，已经编译好，不需要自己编译。
wget http://www.clustal.org/omega/clustalo-1.2.4-Ubuntu-x86_64
mview
mview需要下载源码后进行修改
wget https://github.com/desmid/mview/archive/v1.66.ta...
新型冠状肺炎病毒
病毒正式定名为SARS-CoV-2，导致的疾病正式定名为COVID-19。
病毒定名为2019-nCoV。存个档，以备不时之需。
病毒基因组序列。
WHO的实验室检测方式。
病毒病预防控制所公布的qPCR引物探针如下：
** Target 1（ORF1ab）: **
正向引物（F）：CCCTGTGGGTTTTACACTTAA
反向引物（R）：ACGATTGTGCATCAGCTGA
荧光探针（P）：5’-FAM-CCGTCTGCGGTATGTGGAAAGGTTATGG-BHQ1-3’
** Target 2（N）: **
正向引物（F）：GGGGAACTTCTCCTGCTAGAAT
反向引物（R）：CAGACATTTTGCTCTCAAGCT...
根据转录本编号，提取ensembl的cdna序列
本来想通过解析网页来做这个，但是人家ensembl都提供了下载了，还这么做就觉得太蠢了。。。
hg19下载这个文件
wget ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh37.cdna.all.fa.gz
里面有所有的ensembl的cdna转录本。
下载程序
wget https://github.com/pzweuj/practice/blob/master/python/Ensembl/EnsemblTS.py
修改好程序中hg19的路径。
然后参照说明使用即可，例如，查询ENTS00001转录本
python EnsemblTS.py -i EN...
还是来更新一下clinvar吧
主要是annovar中的clinvar数据库，是20190305的版本，太久没有更新了。
以前也写过一篇clinvar的格式调整，但是貌似已经不适用了。按照以往的方式在ncbi的ftp中找到最新的clinvar并且下载下来然后解压。
wget ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar_20200106.vcf.gz
gunzip clinvar_20200106.vcf.gz
然后按照annovar里的格式，撸了个python脚本进行解析
origin = open("clinvar_20200106.vcf", "r")
annovar = open("hg19_clinvar_20200106.txt...
2019总结
2019年的最后一天。
2019年是我家里蹲大半年后重新回归工作的一年。
工作上，现在的公司给了很多接触不同项目的机会，能学到很多。整个分析流程和代码风格都发生了极大的变化。往小处说就是会每个流程都加传参了，往大了说就是模块化的理念越来越深入。后期因为接手同事的数据库，对mysql也了解了更多一点。家里蹲时常用的R倒是被抛下了，python代码越积越多。
生活中，当然是用游戏来充实。今年玩的大作印象深刻的有_只狼_和_塞尔达传说荒野之息_。只狼是上半年玩的，把故事揉碎到各个物品介绍等细节中确实大大的吸引，2019GOTY实至名归。野吹则是国庆时因为不出远门为了消遣入手了switch玩的，不愧是2017年的GOTY。两个作品都是讲故事的好作品。年底时搬了家，也算是换了一个生活环境。...
共计 204 篇文章，26 页。