python大法好，但是坑爹的版本，坑爹的包管理，折腾死大把大把工作时间。anaconda实现了python2和python3和谐共存，同时conda更方便的进行包管理，支持win，linux，mac多系统平台。

还有一些小小的便利，anaconda已经提供了python科学计算常用包，numpy, pandas, matplotlib等等，python notebook也必然是有的，目前已更名为jupyter notebook，输出文档很方便进行交流展示。

清华大学TUNA镜像网站提供了anaconda安装包，https://mirrors.tuna.tsinghua.edu.cn/anaconda/archive/，不必被下载速度折磨死。

anaconda安装

和python一样，anaconda存在2和3两个版本，目前都更新到4.3.1版，但是，，，还有一个但是，选择安装哪一个版本并没有太大区别，可以方便的通过anaconda环境管理切换不同的python版本。

针对unix like用户，anaconda还有一个非常重要的特性，非root用户也能方便的进行安装(默认安装到~/anaconda2/)，添加~/anaconda2/bin到环境变量$PATH，而后即可享受anaconda。

1
2
3
4
5
6
7
8
9
10

#创建名为python3的环境，python版本为3.4，该命令仅安装python基本模块
conda create --name python3 python=3.4
#创建新环境，同时针对新环境安装anaconda
conda create --name python3 python=3.4 anaconda
#激活python3环境
source activate python3
#关闭环境
source deactivate python3
#删除环境
conda remove --name python3 --all

conda安装包

conda于pip类似，对python包进行管理，还能管理python和conda本身

1
2
3
4
5
6
7
8

#列出当前环境中已安装的包
conda list    
#搜索，安装，升级，删除指定特定的包
conda [search|install|update|remove] package_name
#可以使用-n参数，对指定特定环境中的包进行操作
conda list -n env_name
#升级python或者conda，升级python至当前大版本的最新版本
conda update python

anaconda镜像设置

anaconda主机在国外，conda在线安装包速度慢到怀疑人生，这时绝壁选择使用国内镜像网站，首先将仓库镜像加入conda配置，具体命令为：

1
2
3
4

# 添加Anaconda的TUNA镜像
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
# 设置搜索时显示通道地址
conda config --set show_channel_urls yes

或者直接编辑文件~/.condarc

1
2
3
4
5

ssl_verify: true
channels:
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
  - defaults
show_channel_urls: true

然后，你就可以尽情享受python编程的快乐，不被各种包折磨死

ggplot2实现微生物rank abundance绘图，另发现了reshape2能方便的进行dataframe格式整理，长表格，宽表格相互转换

具体代码实现如下：

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
29
30

library(ggplot2)
library(reshape2)

# otu_table    抽平后的otu table，第一列为OTU编号，最后一列为tax注释
otu_table <- read.table(otu_table, sep="\t", header=TRUE)
otu_table <- otu_table[,2:(length(otu_table)-1)]
for(i in 1:length(otu_table)){    # 计算等级相对丰度
  sample <- otu_table[,i]
  otu_table[,i]<- rev(sort(sample/sum(sample)))
}
rank_abundance <- otu_table
rank_abundance$rank <- 1:dim(rank_abundance)[1]    # 添加等级

# 宽表格整理成长表格
dat <- melt(rank_abundance, id.vars = c('rank'), variable.name='sample')
# 整合样品分组信息
group <- read.table(group, sep='\t', col.names=c('sample', 'group'))
dat <- merge(dat, group, by=c('sample'), all=TRUE)
# 移除0值
dat <- dat[dat$value>0,]

g <- ggplot(dat, aes(rank, log10(dat$value), group=sample)) + 
  geom_step(aes(color=group)) + #ylim(-4,0) + xlim(0,2000) + 
  labs(x="Species Rank", y="Relative Abundance") + 
  # 设置y轴指示刻度，针对不同数据，可能需要调整
  scale_y_continuous(breaks=c(-4,-3,-2,-1), labels=c('1e-4','1e-3','1e-2','1e-1')) + 
  # 调整y轴指示刻度显示方法
  theme(axis.text.y = element_text(angle=90, hjust=0.5, vjust=1)) + 

ggsave("rank_abundance.pdf")

结果展示：

linux 登录

远程登录linux服务器具体原理呢，水平还不够，说不清，先简单粗暴上方法。一般采取ssh的方式远程登录，常用ssh客户端有xshell, putty，本地电脑与服务器间数据传输一般是ftp客户端的方式，可以使用xftp, FileZilla。

登录需要帐号，密码，主机地址(ip或可解析网址)，一般ssh登录端口选择22，ftp登录端口选择21。

basic command

1. ls：列出当前路径（或特定路径）下的文件、文件夹

   -l 先是完整信息，包括文件属主、时间、权限等信息

   -a 列出包括隐藏文件在内的所有问价，linux下隐藏文件以.开头

   -h 文件大小以人类友好的方式显示，直观结果为文件大小以k, m, g等代替bit

   -t 安装文件修改时间进行排列

2. pwd：显示当前所在路径

3. cd：进入文件夹，后接相对路径或全路径

4. mv：更名，移动文件

5. cp：拷贝文件

   -L 拷贝源文件为链接时，追踪拷贝原始文件

   -r 针对文件夹拷贝，递归拷贝文件夹下所有子目录或文件

6. ln：创建文件链接

   -s 创建软链接，相比ln创建链接方式，ln -s创建链接不占用硬盘空间

7. mkdir：创建文件夹

   -p 递归创建多级文件夹，eg：mkdir -p /path/to/a/b/c，其中文件夹a原本不存在，创建a同时创建子目录b子目录c

8. rmdir：删除空文件夹

9. rm：删除文件或文件夹

   -r 递归删除，通常用于删除非空文件夹

   -f 强制删除，不显示提示消息

文件操作

1. cat, zcat：zcat支持打开压缩文件

2. head, tail：显示文件前10（默认）行，后10行

   -n num 指定显示行数

3. more, less

4. paste：横向合并两个文本

5. touch：创建文件

6. cut：

   -f n-m 显示文件第n到m列

   -d 设置分隔符，默认‘\t’

7. join

8. vim：文本编辑器

9. wc

   -l 统计文件行数

打包压缩

1. tar：文档打包，通过调用gzip或bzip2压缩

   -c 打包，创建新的tar文件

   -x 解包，从tar文件中提取文件

   -f tar包文件

   -v 显示当前正在处理的文件

   -z 调用gzip进行压缩或解压缩

   -j 调用bzip2进行压缩或解压缩

   -J 调用xz进行压缩或解压缩

2. gzip：压缩和解压gz后缀文件

   -d 解压gzip压缩文件

   -c 压缩或解压结果输出到STDOUT

   -num [1-9] 数字越小压缩比率越低，速度越快，-1等效于-fast，-9等效于-best

3. unzip：解压zip后缀文件

4. bzip2

权限管理

1. chmod：修改文件读写执行等权限，只能修改属主为自己的文件或文件夹

   4, r代表read权限，2, w代表write权限，1, x代表execution权限

   eg: chmod +x file 添加文件执行权限

   ​ chomd 755 file 修改文件权限为属主读写执行，用户组读写，所有用户读写权限

2. chown：修改文件属主信息，eg：chown -R user:group file

   -R 递归修改子文件、目录属主信息

搜索

1. grep, zgrep：文件内容查找

   -f file 查找pattern在file中

   -v 取反，显示不包含pattern的行

   -c 计数，输出包含pattern的行数

2. find：查找文件，指定文件路径下

3. locate：查找文件，所有硬盘内

4. awk：这个命令可以学，详细教程可参考https://github.com/mylxsw/growing-up/blob/master/doc/%E4%B8%89%E5%8D%81%E5%88%86%E9%92%9F%E5%AD%A6%E4%BC%9AAWK.md

5. sed：这个也可以学，详细教程可参考https://github.com/mylxsw/growing-up/blob/master/doc/%E4%B8%89%E5%8D%81%E5%88%86%E9%92%9F%E5%AD%A6%E4%BC%9ASED.md

系统命令

1. top：显示进程

   -c 显示详细命令

   -u user 显示特定user进程

2. ps

3. kill： kill进程

4. free：查看内存

5. history：查看shell执行历史记录

Linux系统命令可以通过man command或者command --help的方式查看详细帮助信息

github看了业界大神li heng序列解析的脚本readfq，智商着急啊，没看懂，倒是新学了一个新功能，python生成器yield。

readfq.py适用于fasta或者单端fastq文件，但是目前测序主要以双端数据为主，现学现卖，实现了一个paried-end数据解析器，代码如下：

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

#!/usr/bin/env python
# --*-- coding:utf-8 --*--

import os
import sys
import gzip

def parse_pe_fastq(fastq1,fastq2,phred=33):
    from numpy import fromstring,byte
    while True:
        name1 = fastq1.readline().partition(' ')[0]
        name2 = fastq2.readline().partition(' ')[0]
        if not name1 or not name2:
            break
        read1, nothing1, qual1 = fq1.readline()[:-1], fq1.readline(), fq1.readline()[:-1]
        read2, nothing2, qual2 = fq2.readline()[:-1], fq2.readline(), fq2.readline()[:-1]
        qual1 = fromstring(qual1,dtype=byte)-phred
        qual2 = fromstring(qual2,dtype=byte)-phred
        assert name1 == name2, 'fastq1, fastq2 is not paired-end'
        yield name1, read1, read2, qual1, qual2

if __name__ == '__main__':
    fq1 = gzip.open(sys.argv[1],'r')
    fq1 = gzip.open(sys.argv[2],'r')
    for i in parse_pe_fastq(fq1, fq2):
        name, seq1, seq2, qua1, qua2 = i
        # filter, stat, split

通过while循环遍历双端文件，四行为一个单元进行数据读取，每四行作为一个完整单元返回（yeild），进行其他操作（# filter, stat, split）， 一个操作结束后迭代器（parse_pe_fastq）继续解析测序数据，for循环中qua1作为一个numpy.ndarray对象，可进行其他一些q20统计、过滤等操作。

生物狗应该都知道、用过blast的吧，出于代码维护等原因，NCBI在原有核心算法基础上重写了blast，使用分别的blastp, blastn, blastx, tblastn, tblastx代替了原来的blastall，软件包名称也变成blast+（与之前老版本blast区分），以下文中blast都指新版blast+。

blast

ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/下载系统对应版本，linux用户可直接下载编译好的版本，解包急用，简单粗暴，windows用户下载ncbi-blast+，和正常windows软件一样安装(建议安装非系统盘，存放数据库文件，占用系统盘空间)，最后添加环境变量就好。

###blast 输出格式

blast支持txt、xml、ASN.1等多种输出格式，具体可通过参数-outfmt int参数进行设置（int取值0-11 ），其中以xml（-outfmt 5 ）和txt（-outfmt 6 ）使用较多。

xml输出信息全面，可能不便于查看，txt格式便于查看，而仅显示有限的匹配信息，但是blast支持自定义输出信息，可通过-outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore"实现，具体输出信息参考blastp -help信息。

###blast task 参数设置

一直很好奇为啥网页版 blast 中的 blastn 有megablast, discontiguous megablast, blastn三种选项，而在单机版中找不到对应可执行程序，看NCBI介绍也有说 megablast 对比速度上更快一些，一直想要自己动手试试，找不到执行程序。某天瞎逛各种论坛的时候，偶然发现了还有-task这么一个参数选择不同的比对模式，就先记录了解下，准备后续测试。

• blastp，-task参数有blastp , blastp-short和deltablast三种选择
• blastn，-task参数可选blastn, blastn-short, dc-megablast, megablast, rmblastn，这里blastn软件默认就采用megablast进行同源搜索
• blastx，-task参数可选blastx, blastx-fast，测试结果可以参看 genedock 博客，https://www.genedock.com/blog/2017/02/22/20170222_blastx-survey/

###其他比对工具

1. Diamond, 类blast程序吧，运行时间、速度上甩blast+ n条街，空间换时间的策略吧，占用内存、CPU也是够够的了，不过最新版本加入了-b参数限制内存、虚拟磁盘空间的占用，同时新增--sensitive和--more-sensitive参数比对结果的灵敏性
2. Blat, UCSC出品，没仔细研究，某些时候比 blast 好用

生物背景从事生物信息分析工作，本科时候对统计、计算机之类的有些兴趣，觉着“生物信息”很高大上的样子，什毛都不懂就如了生信的坑，这一路下来，知道越多越觉得需要了解、懂的东西也越多，经常有“我擦勒，要学这么多东西啊” ， linux, R, python, perl…… 微博上看了各种从入门到放弃系列，觉得很有意思，这里也不自知的写个从入门到放弃吧，生物信息是一个多学科的交叉学科，遗传、生理、计算机、统计基础好像都须知道一些，“XXX，从入门到放弃”，我也不知道是啥从入门到放弃。

很早就有利用博客整理、交流学习思路的想法，懒癌穷癌很成功的阻止了我，看到很多利用github的教程，简单粗暴，终于鼓起勇气试试。

讲真，统计不是很懂，遗传不知道怎么说起，那就说说生信分析的linux、R、python这些咯。

Linux

分析主要工作环境，命令行操作，分析工具众多。

基础命令：cd，pwd，ls，mkdir，rm，touch，cp，mv，echo，wc，ln

文件查看：head，tail，more，less，cat (zcat)

打包压缩：tar，gzip，zip，bzip2

文件处理：vim，awk，sed，grep (zgrep)，join

其他几个有用的命令：find，top，ps，kill，free

Linux系统命令可以通过man command或者command --help的方式查看帮助信息

其他几个须知的点：重定向>，文本追加>>，管道|，后台运行&和nohup，标准输入、标准输出和标准错误

R

两个功能，分析和绘图。分析这个主要是利用一些现成R package，例如转录组分析中DESeq、edgeR这些；绘图，当时还是ggplot2，教程很多，我这个R战五渣没有发言权。学习R，当然要有RStudio，R书籍什么的， R语言实战、ggplot2数据分析与图形艺术都可以看下。

python

主要是文件整理吧，不同格式转换，提取有效信息。python语法简单，基本上很快都能写出实现自己功能的代码，有很多第三方modules帮你实现想要的功能，另外学习的人多啊，不怕找不到人问问题。

编程能力这个还是靠多练习了，一些知识点就只能自己看书看教程背，现在网上也很多python相关的资料教程，边看边练习，入门很快的。记得，我当时看了很多七七八八的教程，然而少练习，折腾挺长时间，廖雪峰的python教程 (2.7)，笨方法学python ，github上也是有很多python学习笔记或者英文书籍翻译。

强烈推荐ipython，可以理解为python里的Rstudio吧，可以生成富文本(filename.ipynb)文件，很方便用于脚本展示、交流。

Anaconda ，一个python科学计算环境，包括了常用的numpy、pandas、matplotlib、ipython等modules，管理python工作环境，实现python2, python3和谐共存，这个也是刚开始了解，熟悉的朋友可以交流交流。

perl

和python类似吧，文件整理。这个这个，只会打印hello, world，算是不会了，网上也很多博客教程，可以参考学习。

其他技能

MarkDown

语法简单，结果展示酷炫，低门槛装逼利器。

Typore：免费markdown编辑器里面最好的，界面简洁，支持多种格式（pdf，html，docx等等）导出。

Docker

可以掌握的技能。关注docker，这是因为同事提到工作流程迁移系统环境兼容的问题，docker完全免去了这些烦恼。无奈，工作系统版本太低，无法安装，还没有试过。