Python处理临时文件

tempfile这个模块主要是用来创建临时文件和目录,用完后会自动删除,省的你自己去创建一个文件、使用这个文件、再删除这个过程了。其中比较常用的是TemporaryFile和NamedTemporaryFile。

import tempfile
import os
 
tmp1 = tempfile.TemporaryFile()   # 创建一个临时文件
type(tmp1)   # tempfile._TemporaryFileWrapper
tmp1.close()    # 关闭后,文件立即删除
 
tmp2 = tempfile.NamedTemporaryFile(delete=False) # 创建有名字的临时文件,且关闭后不立即删除
tmp2.name    # 'C:\Users\chenw\AppData\Local\Temp\tmpj9whdpr8'
os.path.exists(tmp2.name)   # True
tmp2.write("This is a tmp file!".encode())    # 不接收字符串(str),只能是字节流(bytes)
 
tmp2.seek(0)    # 调用read()函数前,需要调用seek()函数,指定读取指针的位置,参数是偏移量,0表示从头开始
tmp2.read()    # b'This is a tmp file!'   read()函数的返回值类型依然是bytes,而不是str
tmp2.close()    # 关闭后,文件不立即删除

Python发送邮件

Python提供了标准库smtplib用于发送邮件,整个过程非常简单,不足的地方是密码是明码。

import smtplib
from email.mime.text import MIMEText
 
#发件人
mail_from = "BioChen.com<admin@biochen.com>"
mail_host = "smtp.exmail.qq.com"
mail_usr = "admin@biochen.com"
mail_pwd = "pass_word"
 
#收件人
mail_to_addr = "pm.chenwen@qq.com"
 
#登录服务器
server = smtplib.SMTP(mail_host)
server.login(mail_usr,mail_pwd)
 
#邮件内容
content = 'this email is a test for python!'
msg = MIMEText(content,_subtype='plain')
msg['Subject'] = 'this email is a test for python!'
msg['From'] = mail_from
msg['To'] = mail_to_addr
 
#发送邮件
server.sendmail(mail_from, mail_to_addr, msg.as_string())
 
#关闭服务器
server.close()

网易和腾讯邮箱的SMTP服务器地址:

服务商 发件服务器 端口
网易163邮箱 smtp.163.com 25
网易企业邮 smtp.ym.163.com 25
腾讯QQ邮箱 smtp.qq.com 25
腾讯企业邮 smtp.exmail.qq.com 25

Ubuntu 18.04配置固定IP

Ubuntu 18.04不再使用ifupdown配置网络,而改用netplan。在/etc/network/interfaces配置固定IP是无效的,重启网络的命令services network restrart或/etc/init.d/networking restart也是无效的。

1. 使用ifupdown配置网络
如果要使用之前的方式配置网络,需要重新安装ifupdown:

sudo apt install ifupdown

修改配置文件/etc/network/interfaces:

sudo vim /etc/network/interfaces

配置文件修改如下:

iface ens160 inet static
address 210.72.92.25
gateway 210.72.92.254
netmask 255.255.255.0
dns-nameservers 8.8.8.8

重启网络服务使配置生效

sudo service network restrart

2. 使用netplan配置网络
Ubuntu 18.04使用netplan配置网络,其配置文件是yaml格式的。安装好Ubuntu 18.04之后,在/etc/netplan/目录下默认的配置文件名是50-cloud-init.yaml,我们通过VIM修改它:

sudo vim /etc/netplan/50-cloud-init.yaml

配置文件修改如下:

network:
    ethernets:
        ens160:
            addresses:
                - 210.72.92.28/24 # IP及掩码
            gateway4: 210.72.92.254 # 网关
            nameservers:
                addresses:
                    - 8.8.8.8 # DNS
    version: 2

无论是ifupdown还是netplan,配置的思路都是一致的,在配置文件里面按照规则填入IP、掩码、网关、DNS等信息。注意yaml是层次结构,需要缩进,冒号(:)表示字典,连字符(-)表示列表。

重启网络服务使配置生效:

sudo netplan apply

搭建远程Jupyter Notebook服务器

Jupyter Notebook是数据分析的一大神器,特别是可视化和分享功能强大。数据量大的时候,我的笔记本电脑有点吃力,于是想着在服务器上搭建Jupyter Notebook服务器。

1. 安装Jupyter Notebook库

$ pip3 install Jupyter

2. 生成Jupyter Notebook配置文件

$ jupyter notebook --generate-config
Writing default config to: /home/you/.jupyter/jupyter_notebook_config.py

3. 设置Jupyter Notebook密码

$ jupyter notebook password
Enter password: *******
Verify password: *******
[NotebookPasswordApp] Wrote hashed password to /home/chenwen/.jupyter/jupyter_notebook_config.json

4. 修改jupyter_notebook_config.py

c.NotebookApp.ip = '*'               #所有绑定服务器的IP都能访问
c.NotebookApp.port = 9999            #将端口设置为自己喜欢的吧,默认是8888
c.NotebookApp.open_browser = False   #我们并不想在服务器上直接打开Jupyter Notebook,所以设置成False

5. 启动Jupyter Notebook服务器,并让其在后台运行

$ jupyter notebook &

6. 在自己的电脑上打开浏览器访问Jupyter Notebook服务器

输入密码后,就能像本地Jupyter Notebook一样使用啦。

更多细节请参考Jupyter Notebook官方文档:Running a notebook server

Google翻译API

Github上面其实有Google翻译的API,由于我国社会主义初级阶段的基本国情,直接拿来用好像行不通!

import requests
from bs4 import BeautifulSoup
 
def getHTMLText(url):
    try:
        r = requests.get(url, timeout=30)
        r.raise_for_status()
        return r.text
    except:
        print("Get HTML Text Failed!")
        return 0
 
def google_translate(to_translate, from_language="en", to_language="ch-CN"):
    #根据参数生产提交的网址
    base_url = "https://translate.google.cn/m?hl={}&sl={}&ie=UTF-8&q={}"
    url = base_url.format(to_language, from_language, to_translate)
 
    #获取网页
    html = getHTMLText(url)
    if html:
        soup = BeautifulSoup(html, "html.parser")
 
    #解析网页得到翻译结果
    try:
        result = soup.find_all("div", {"class":"t0"})[0].text
    except:
        print("Translation Failed!")
        result = ""
 
    return result
 
print(google_translate("Hello World!"))
 
#你好,世界

代码好像超级简单哈,有空的话我想把它做成一个库。有github的给个星星呗,鼓励一下!
https://github.com/wen-chen/PyTransl

Python图像库Pillow

PIL(Python Imaging Library Python图像库)是一个强大的图像处理标准库,功能强大却又简单易用。现在的名字叫做Pillow。

安装Pillow

pip3 install pillow

如果你使用的是Anaconda开发环境,自带Pillow库。值得注意的是,导入Pillow库时,库的名称是PIL。下面我们将用Pillow库做一些有意思的事情。

1. 图像手绘效果

from PIL import Image
import numpy as np
 
a = np.asarray(Image.open('./picture.jpg').convert('L')).astype('float')
 
depth = 10. 					#(0-100)
grad = np.gradient(a)				#取图像灰度的梯度值
grad_x, grad_y = grad 				#分别取横纵图像梯度值
grad_x = grad_x*depth/100.
grad_y = grad_y*depth/100.
A = np.sqrt(grad_x**2 + grad_y**2 + 1.)
uni_x = grad_x/A
uni_y = grad_y/A
uni_z = 1./A
 
vec_el = np.pi/2.2 			        #光源的俯视角度,弧度值
vec_az = np.pi/4. 			        #光源的方位角度,弧度值
dx = np.cos(vec_el)*np.cos(vec_az) 	        #光源对x 轴的影响
dy = np.cos(vec_el)*np.sin(vec_az) 	        #光源对y 轴的影响
dz = np.sin(vec_el) 			        #光源对z 轴的影响
 
b = 255*(dx*uni_x + dy*uni_y + dz*uni_z) 	#光源归一化
b = b.clip(0,255)
 
im = Image.fromarray(b.astype('uint8')) 	#重构图像
im.save('./picture2.jpg')

图像效果如下:

2.生成数字验证码

from PIL import Image, ImageDraw, ImageFont, ImageFilter
import random
 
#随机字母:
def rndChar():
    return chr(random.randint(48, 57))
 
#随机颜色1:
def rndColor():
    return (random.randint(64, 255), random.randint(64, 255), random.randint(64, 255))
 
#随机颜色2:
def rndColor2():
    return (random.randint(32, 127), random.randint(32, 127), random.randint(32, 127))
 
#240 x 60:
width = 60 * 4
height = 60
image = Image.new('RGB', (width, height), (255, 255, 255))
#创建Font对象:
font = ImageFont.truetype('ariblk.ttf', 40)
#创建Draw对象:
draw = ImageDraw.Draw(image)
#填充每个像素:
for x in range(width):
    for y in range(height):
        draw.point((x, y), fill=rndColor())
#输出文字:
for t in range(4):
    draw.text((60 * t + 10, 10), rndChar(), font=font, fill=rndColor2())
#模糊:
image = image.filter(ImageFilter.BLUR)
image.save('code.jpg', 'jpeg')

得到的数字验证码如下图:

windows中安装Tensorflow

TensorFlow 是一个用于人工智能的开源神器,我很想学习它。TensorFlow目前在windows平台上面只支持Python 3.5,我的电脑上面装的是Anaconda3-4.3.0.1,带的Python最新版本3.6。我不想将Python版本降级,于是就在虚拟环境中安装TensorFlow吧。

1. 创建TensorFlow虚拟环境
Win+R,运行”CMD”,输入以下命令

conda create --name tensorflow python=3.5

这样我就创建了一个名叫tensorflow的虚拟环境,当然你可以用任何你喜欢的名字代替。我还指定了Python版本是3.5。

2. 激活虚拟环境

activate tensorflow

如果想取消激活,使用deactivate命令即可

deactivate tensorflow

3. 安装TensorFlow

conda install tensorflow

Anaconda会自动安装TensorFlow和依赖关系

4. 安装Spyder和Jupyter Notebook
我非常喜欢用Spyder和Jupyter Notebook写Python代码。

conda install spyder
conda install jupyter notebook

这样我们的Windows TensorFlow开发环境就创建好了!

Win+R,运行”CMD”,

activate tensorflow
spyder

打开的是Spyder(Python3.5),而不是真实环境下的3.6。在交互窗口运行import tensorflow也没有报错,说明我们环境搭建成功。

更多Anaconda的用法参看我之前写得一篇博客:Python科学计算工具箱Anaconda

佛祖的灯芯

《大话西游》里有一段:

紫霞来到了驿站的外面,在院落里趴着一只黑狗,正在盯着紫霞,一个老者正在来回的奔波着,不断的在院子里取火,然后又跑进屋里去点灯。
“各位客官,请稍等一下,马上就好了。”
灯里面只有灯油,没有灯芯,当然没办法点着了,紫霞走了过去:“老伯啊,你这样是点不着的,你只有油,没有灯芯哎。”
“我知道,可是灯芯不见了,有什么办法呢?”
“再买一根喽。”
紫霞和二郎神化作的老头说了几句,老头忽然直起身来,看着紫霞,声音骤然一沉。
“你说得对,可是如来佛祖的那盏日月神灯的灯芯不见了,到哪里去买?”

没有灯芯点不着灯,那人心呢?

PhastCons、PhyloP和保守性

全基因组多序列比对的主要目的是揭示基因组中进化上保守的和非保守的区段以及它们的分布。但是,在碱基层次观察多个全基因组的多序列比对结果和序列保守性非常不方便;其次,许多功能区段的保守性介乎高度保守和完全不保守之间;再者,保守性要能予以方便的显示。上述三点使UCSC基因组浏览器使用两个软件PhastCons和PhyloP把由MULTIZ产生的多序列比对结果转换成两种单一的保守性计分和显示,这两个方法都基于已知的种系树结构和利用一个称作phylo-HMM的隐马尔可夫模型种系分析方法。

PhyloP只考虑比对的当前列而PhastCons同时也考虑比对列的相邻列,这使得PhastCons对于保守区段的出现更敏感而PhyloP对于保守区段的界定更有效。PhyloP和PhastCons之间的另一个区别是,PhyloP能够识别加速进化和保守进化的位点,它们分别产生正的和负的计分,而PhastCons的计分总是一个0至1之间的正值。

PhyloP和PhastCons多物种比对文件可以从UCSC基因组浏览器的网站上下载,他们还有提供了bigWigAverageOverBed这个工具。

用法非常简单,只是输入文件in.bw、in.bed和输出out.tab必须是这个名字。

使用MCScanX分析基因组共线性区块

MCScanX采用改进了的MCScan算法,分析基因组内或者基因组间的共线性区块。它利用两个物种蛋白质blastp比对结果,再结合这些蛋白质基因在基因组中的位置,得到两个物种基因组的共线性区块。如果是分析基因组内的共线性区块,物种内蛋白质自己比对自己就好了。

一、MCScanX安装
下载MCScanX的安装包MCScanX.zip,可以在linux系统和Mac OS进行编译:

$ unzip MCscanX.zip
$ cd MCScanX
$ make

如果是在64位系统进行编译,msa.h、dissect_multiple_alignment.h、detect_collinear_tandem_arrays.h这三个源文件的开头需要添加#nclude
MCScanX包括MCScanX、MCScanX_h、duplicate_gene_classifier三个主程序,位于主文件夹中;还有12个下游分析程序位于downstream_analyses文件夹中。MCScanX 检测共线性区域,并比对到参考染色体上。MCScanX_h 和MCScanX类似,只不过输入文件是成对的用tab分隔的同源基因。下游分析程序主要用于可视化输出。

二、输入文件的准备
我们需要两个物种的蛋白信息,还有对应蛋白在基因组中的位置。推荐从Ensembl下载蛋白质信息,它的注释信息非常全。如下是Homo_sapiens.GRCh38.pep.all.fa(下载于Ensembl,人的全部蛋白质信息)的一条蛋白质的信息:
>ENSP00000452494.1 pep chromosome:GRCh38:14:22449113:22449125:1 gene:ENSG00000228985.1 transcript:ENST00000448914.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:TRDD3 description:T cell receptor delta diversity 3 [Source:HGNC Symbol;Acc:HGNC:12256]
该蛋白质名ENSP00000452494.1,位于14号染色体上,起始位置22449113,终止位置22449125。

blast建库:

$ makeblastdb -in a.fa -dbtype prot -parse_seqids -out adb

blast比对:

$ blastp -query b.fa -db adb -out xyz.blast -evalue 1e-10 -num_threads 16 -outfmt 6 -num_alignments 5

按照软件说明书的要求,-evalue 1e-10是将e value卡在1e-10,-num_alignments 5是取最好的5个比对结果,-outfmt 6是输出格式为tab分隔的比对结果。

蛋白质基因在基因组上的位置文件为gff文件。这里的gff文件不是平常用到的那种,需要修改一下。文件是一个tab分隔的,只有4列:
sp# gene starting_position ending_position
sp是两个字母,后面的#号表示染色体编号,这个文件自己生成一下好了。

三、使用MCScanX分析基因组共线性区块
假设我们的blast文件为xyz.blast,gff文件为xyz.gff,且都放在一个文件夹下,我们可以简单运行一下命令:

$ MCScanX xyz

程序输出的结果有两个xyz.collinearity文件和xyz.html目录。
xyz.collinearity,成对的共线性区域,如下图所示:

xyz.html目录里面有很多小文件,文件名称是参考基因组染色体编号。第一列是每个基因位点的复制深度,第二列是基因参考染色体,红色部分是串联基因,后面黄色部分的内容是比对上的共线性区域,只有哪些比对上的基因会被展现出来。html文件可以用网页浏览器打开。如下图所示:

下面是一些更详细的参数:

四、下游分析程序
下游分析程序很多,主要是用来出图的,就不一一介绍了。
dot_plotter.java画散点图

$ java dot_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

circle_plotter.java画圆圈图

$ java circle_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

输入除了MCScanX的结果和gff文件,还需要一个ctl文件,这个文件自己按照MCScanX的说明去写即可,告诉程序图要画成什么样子。

更多用法参见官方的说明书。
http://chibba.pgml.uga.edu/mcscan2/documentation/manual.pdf