mash 是一个可快速计算两个（或多个）基因组之间相似度的软件，可输入基因组序列或测序结果（fastq文件），使用计算机科学中一个经典的方法MinHash，一般用来判断两个网页内容是否雷同，以决定是否从搜索结果中去掉某些类似的结果。

使用示例：

mash dist -p 4 -l 、/path/to/ref.fa query_fa.txt > query_mashed.txt

用于计算所有query序列和ref序列的相似度，而相似度和物种间的亲缘关系是相关的，其中<query_fa.txt> 是一个记录了所有query序列的文件

ABI测序仪下机的结果包含.ab1和.seq文件，通常可以导入BioEdit这样的图形化软件检查峰形和运行blast，而如果想要实现自动化的检测杂合位点并做比对，最终得到变异位点，就需要花些时间研究下。

主要的问题是，ab1文件是二进制的，关于文件格式的文档稀少，唯一靠谱的是ABI的官方文档 仔细看过这个文档你会发现，看完之后作为新手的你尽管对于这个格式有所了解，知道它分为Header、Directory以及Data，以及有很多的tag。但知道这些对于想要从ab1提取trace value并进行heterogeneous base calling是远远不够的。所好的是，已经有人有过类似需求并作了尝试。

• abifpy 实现的一个解析ab1文件的模块，现已并入BioPython。优点是兼容Python2和Python3，部分meta data已经给你提取出来了，而如果想要更多，可以使用get_data(‘youkey‘)获取。

• 如果只需要取出每个通道的trace value，有了官方文档和abifpy就够了，如果还想学习更多，可以读一读SeqTrace源码，从中可以学到怎样组织一个项目的代码，以及怎么写测试。

• 前面两个都研究一番之后，还是不能满足你的好奇心？可以把io lib这个库研究下，学习使用它的API.

翻了一下以前的文章发现，最近的上一篇文章是2013年11月的，那时候还在南京的学校，从找工作开始，然后面试，来北京，上班之后就更没有留出时间静下来写东西了。后来换了Mac，还特意把blog的资料作了备份，想着不久之后，重搭环境另起炉灶的，没想一下子拖到不久前才在新电脑上配好环境。

上班的这一年多，改变了很多。工作上从项目执行到流程负责人，再到信息分析主管而后技术副总监，见证了公司规模的扩张和变迁，也目睹了身边同事的来来去去，“送走”了三波运营岗同事，现在接触的大都是第四波的人。公司在成长，自己的工作也在经历变化，以前只需要自己干活，现在要考量人和事，把相应的事交给合适的人去做，然后跟踪进度，培养新人。

除了工作，生活上也有一些变化。刚到北京住的没有独立卫生间的小院，后来的院子有单独卫生间，尽管很少，还是方便多了。随公司搬到天津之后，住的更好。经过两段感情经历的洗礼，也成长了很多，对自己也有了更多地认识。写东西是一个特别好的总结的方式，通过它可以和自己对话，从而自省；和别人交流，从而结交同好。所以，以后我会坚持下去的。

I found this when hanging around online. To be honest, I am dying to learn web crawling long time ago, can not wait to play with this.

modules involved:

• urllib2
• re
• BeautifulSoup Need to install manually

This post is very easy to follow with basic familarity of Python, the difficulty of web crawling is lots of, if not all, websites follow HTML specification exactly, there are often errors here and there if you try to parse useful bits of infomation out there, and here comes BeautifulSoup which handle most of these for us.

The second section discuss sending email with Python. One thing need to keep in mind is check your internet connection first if your code refuse to work :-)- I have complete the spec of No.5, and below are some commands I find useful:

• execfile('script.py') #load script into console
• python -i script.py #execute script then come to Python console, ‘i’ short for ‘interactive’
• sometimes you need to insert a time stamp into you file, and time.ctime() is probably the simplest solution
  source code is here

前几天接到一个公司的电话面试，当天晚上HR发来一道题，限三日内完成并发回。我看了下，正好用到的知识自己都熟悉，花了几个钟头就给发回去了 :-)

原题大意这样 写一个程序，分别计算一个目录里所有fasta格式的核酸序列的GC含量，然后按照GC含量由低到高显示，并把结果写入一个文件
显示序列信息的格式为：
文件名 序列名 GC含量
hem.fasta hem1 24.3%

这道题挺简单的，细看之下主要有这么几点：

• 计算序列GC含量，可以用一个函数实现
• 列出目录下的所有文件，然后判断他们是不是fasta格式文件。前者用os模块的listdir就可，后者就是要看文件的首行是不是形如“>at1g23490”开头，即是否“>”开头，后面紧跟字符，正则就行了。
• 结果排序。这个应该是这道简单笔试题最难的部分了。正好之前从一个博客读到解决排序的办法，就依样做了。把每条序列的结果存到一个列表里，然后对列表排序。

这是我的代码 信息分析面试

EOF

###git
在github上参与别人的项目或则多人共同开发一个项目的话，fork一个项目的仓库过来，自己做完提交pull request之后，就等着作者review过之后进行merge了。但是之后作者再进行新的改动，我的分支上还是原样，那么怎样让把我的分支上的与最新的仓库一直呢?

1
2
3
4
5

#assume you originally cloned from your forked repo,that remote will be called "origin", 
#then you need to add firstguy repo as another remote
git remote add firstguy git://github.com/firstguy/repo.git
#After that is all set up, you should indeed be able to
git pull firstguy mastergit push origin

###ubuntu apt
apt search

1
2
3

# package search
apt-cache search 
apt-cache show #------(package 获取包的相关信息，如说明、大小、版本等)

apt install&remove

1
2
3
4
5
6
7

#package install
sudo apt-get install 
sudo apt-get install # -----(package - - reinstall 重新安装包)
sudo apt-get -f install # -----(强制安装?#"-f = --fix-missing"当是修复安装吧...)
sudo apt-get remove #-----(package 删除包)
sudo apt-get remove - - purge # ------(package 删除包，包括删除配置文件等)
sudo apt-get autoremove --purge # ----(package 删除包及其依赖的软件包+配置文件等

apt update and more

1
2
3
4
5
6
7
8
9
10
11

#update source
sudo apt-get update 
sudo apt-get upgrade #------更新已安装的包
sudo apt-get dist-upgrade # ---------升级系统
sudo apt-get dselect-upgrade #------使用 dselect 升级
apt-cache depends #-------(package 了解使用依赖)
apt-cache rdepends # ------(package 了解某个具体的依赖?#当是查看该包被哪些包依赖吧...)
sudo apt-get build-dep # ------(package 安装相关的编译环境)
apt-get source #------(package 下载该包的源代码)
sudo apt-get clean && sudo apt-get autoclean # --------清理下载文件的存档 && 只清理过时的包
sudo apt-get check #-------检查是否有损坏的依赖

nohup: short for “no hang up”

• usage: nohup command &
• default output: nohup.out which contains std out and std error, I refquently use this: nohup sslocal > log &

EOF

Here is the issue: When I am trying to help someone to reinstall wondows7 on a TOSHIBA laptop, I can not get into the GUI interface, no matter with windows PE, USB disk created by universal USB installer or windows 7 DVD download tool which is released by microsoft.

After all these failed, I turn to ubuntu. I plug in my USB disk with ubuntu 12.04 inside, press F12 when booting, and I finally can have my desktop where I can browe all files on the hard drive. I had to erase the system partition when the attempt to check and repair the possible error in system partition failed me.

Then all things go as expected, so you see, that is one reason to learn linux :-) -

EOF

从rna-seqblog看到这篇paper,他们对流行的寻找差异基因的方法多了评测，压根没提到DEGSeq,我看到在bioconductor上面，DEGSeq的下载量只有DESeq、edgeR的一半，所以人家说‘流行的’就把DEGSeq排除在外了。

不过这里要说的主角是另一本书，在翻看前述论文时看到的，打来不禁有相见恨晚的感叹，原书还没看完，这里把里面的一些要点记录下来。
###microarray和sequencing
RNA-Seq相比microarray的优点已经被谈得够多了，这里作者同样讲了RNA-Seq的不足：

• RNA-Seq数据有GC偏好
• mapping的不确定对基因组不同区域的影响是不一样的
• RNA-Seq比较不同处理下的基因表达水平有偏好性。RNA-Seq的统计方法对高counts的检测能力要比低counts强，也就是说RNA-Seq会倾向于将长度更长的基因算作差异（DE）基因。
• RNA-Seq的高灵敏性带来的问题是，提取、富集RNA，RNA片段化和转化成cDNA的过程中都有更大的带来偏差的可能。如小RNA提取的方法强烈影响得到的序列集合。
• RNA-Seq数据更复杂，数据分析可能是瓶颈。

###mapping
主要有两类方法：基于hash-tables;基于Burrows Wheeler Transform

• hash-tables优点是可以找出用户指定的所有structural variants,确定是内存要求高。
• BWT相比前者使用内存更经济，缺点是，BWT确认错配的方法涉及把每个read的大量可能组合进行比对，二这需要相当大的计算量。很多实现方法对此的处理是，不穷尽所有错配比对，这也就意味着一些正确的比对丢失了。

###RNA-Seq实验设计

• 生物学重复重要性：如果我们比较一个癌症组织和一个正常组织，那么我们下的任何一个结论都不能推广到其癌组织群和正常组织群，因为分析过程没有考虑癌组织与癌组子和正常组织与正常组织间潜在的巨大差异。我们的结论仅适用于该研究的样本。
• 上样的时候要考虑到不同lane，不同flow cell之间存在的系统差异
• reads长度的确定。更长的reads花费多，但map到trancriptome会更多，对外显子-外显子间连接区分能力更强，发现更多的遗传变异，但是DE的发现来自reads密度，而不是序列覆盖度。总之，是用更长的reads来换取‘mappability’，还是用更深的测序深度换取更多的reads密度呢？作者给的答案是，投资与‘mappability’所获通常比较小。
• paried-end reads: PE-reads可以有更多的reads map到转录组或基因组上，有助于发现结构变异，然而PE-reads带来的更多花费并不能增加检测DE的强度。

EOF

Our DGE project report provide a differently expressing genes list which is infered using DEGseq. When I was trying to follow the analysis pipline, I found that DEGseq has a so simple document that I cannot even figure out how to apply the function in my work, then I turn to DESeq, which provide a very good document with explanation of options and typical cases. Now I have to DG(Different Gene)lists: one from NOvogene, another is mine(generated by DESeq). Some genes are the same in them, while some not. My goal is to identify the same and plot a venn diagram. Here is the outline:

• idnetify the same genes, calculate the number of same and different genes.
• draw a venn diagram in R to display above result.

DEGseq list has 2k+ genes, and my DESeq list has 1230 genes, obviously I need a quick script to compare each gene in DESeq and remember the number of genes in both.

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
31
32
33
34
35
36
37
38
39
40
41

import sys



dfgnlist1 = sys.argv[1]

dfgnlist2 = sys.argv[2]

with open(dfgnlist2, 'r') as deg:

  dfgn2 = deg.readlines()

with open(dfgnlist1,'r') as de:

  dfgn1 =de.readlines()
i = 0

j = 0

for g in dfgn1:

  if g in dfgn2:

    with open("samegenes.txt",'a') as samegenes:

      samegenes.write(g)

      i += 1

  else:

    with open("dfgenes.txt",'a') as dfgenes:

      dfgenes.write(g)

      j += 1
print '========'

print i, 'same genes'

print j, 'different genes'

用tophat做完mapping之后，接下来就是count reads，HTSeq的文档很详细，还有一个use case，除了这个case，其他的全部跟着走了一遍。接下来使用htseq-qa和htseq-count的时候，遇到了问题，记录在这里：

###htseq-qa
安装了HTSeq后，htseq-qa的使用也很简单：

• htseq-qa [options] read_file #assumes htseq is in your PATH
• python -m HTSeq.scripts.qa [options] read_file #if htseq is not in your PATH

检查clean data的时候很顺利，在read_file是Col_0的raw data时，遇到了问题，原因是raw data里面有的reads quality string信息不完整，尝试过删除该reads或者把reads quality数据补全，都不成。

###htseq-count
htseq-count的使用与qa类似，如果htseq不再你的PATH，可以这样：python -m HTSeq.scripts.count [options] <sam_file> <gff_file>。我遇到的问题有两个：

1. tophat输出的是BAM格式的结果，而htseq-count需要SAM格式的比对结果，首先要转换，这个用samtools就行了：samtools view col.bam >> col.sam samtools默认的输入就是BAM格式，这里用重定向将屏幕输出写入col.sam
2. 然后运行htseq-count报错：”ValueError: The attribute string seems to contain mismatched quotes”.搜了下，报错的行里面’gene_name’id含有多余的分号，手工修改了几次，还是报同样的错。这时候需要写一个脚本了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

import re
import sys

input_file = sys.argv[1]
output_file = sys.argv[2]
j = 0
with open(input_file, 'r') as f:
  for line in f.readlines():
    gene_name = re.search(r'gene_name\s(\S+)',line)
    gene_id = gene_name.group(1)
    match = re.search(r'(\w+);(\w+)', gene_id)
    if match:
      newid = ''.join(['"', match.group(1), match.group(2), '"', ';'])
      newline = line.replace(gene_id, newid)
    else:
      newline = line

    with open(output_file, 'a') as newf:
      newf.write(newline)
      j += 1

print 'totaly', j, 'number of lines modified!'

得到正确格式的gtf之后，就可以count了：python -m HTSeq.scripts.count -m intersection-nonempty -s no sam_file gtf_file >> samplecount 需要注意的事这里’-s’选项默认是’yes’,所以如果样品建库不是连特异性的话，一定要加-s no,要不有一半的reads就浪费了。另外就是把count结果写入文件，以备之后分析。