水一篇。我喜欢这种类似回文串的日期。 有三个多月没更新博客了,马上满血复活,要填的坑还是比较多,希望能早日填完。

坑之列表

近期

  • animation package bug
  • Freshman21 theme bug
  • 欠的blog文章 (10+)
  • 某bot

远期

  • 某R包编写

久远

  • 某专题

其余部分包括Git的使用,可重复性研究(Rmarkdown),以及一些简单的生物信息学分析(Blast,DNA mutation,RNA expression),内容不连贯,我就不一一罗列了。

我当时有个笔记草稿,现已放到了gist上,欢迎参观,但应该没人能完全看懂。。。

2017年更新

最新知识,比较RNA表达量的分布,需要知道最少用几个样本,有了样本才能知道表达量的分布,那么究竟表达量分布长什么样子,需要多少个样本?


在R语言部分,老师推荐了这本书An Introduction to Statistical Learning with Applications in R

还稍微介绍了一点编译相关的内容。

inHybrid Hybrid compiled
Python Java C
R Perl C++
Bash    

接下来就是简单的实战,我发现没有编成经验的学生最容易弄错的一段代码是随机数,他们一般都会问,为什么自己的结果同老师演示的结果不同。

x=100
myfunction1<-function(x){
  v=rnorm(x)
  mean(v)
}
myfunction1(x)

全局变量

在课堂中老师忘了如何将函数中的变量定义为全局变量,其实这个功能我在R中基本就没有用到过。不过我还是查到了方法,用双箭头<<-来定义全局变量。

x=100

myfunction=function(x){
  v<<-rnorm(x)
  mean(v)
}
myfunction(x)
v

杀掉Rstudio进程

有时由于程序写错了,进入了死循环,并且内存占用越来越大,最简单的处理方式是在终端中杀掉Rstudio进程(虽然不推荐)。

killall rstudio

R的数据结构

  • vector
  • list
  • data.frame (还有一种data.table)
  • matrix
  • character

课上练习

words=c("Hello","and hi","World")
l=list()
for(i in 1:length(words)){
  if(i==3 && (words[i]=="World" || words[i]=="world")){
    l[[i]]="America!"
  }else{
    l[[i]]=words[i]
  }
}
for(i in l){
  cat(i," ")
}

实用网站和关键函数介绍

正则表达式

没教\w,老师的观念是用\S,\s,\d,\d+来推导一切表达式。


fh=file("/media/sf_CB2/dm/dm.faa",open="r")
while (length (line=readLines(fh,warn=F,n=1))>0){
  pattern="^\\>gi\\|(\\d+)\\|ref\\|(\\S+)\\|"
  #m=regexec(pattern,line,perl=T)
  #m=mymatch(pattern,line)
  #if(m[[1]][1]==-1){
  #  next
  #}
  #if(!(ifmatched(m))){
  #  next
  #}
  m=mymatch(commapattern,line)
  if(!ifmatched(m)){
    my=mymatch(pattern,line)
  }
  if(!(ifmatched(m))){
    next
  }
  cat(getstring(line,m,2))
  cat("\n")
  #v=regmatches(line,m)
  #vec=v[[1]]
  #cat(vec[2],vec[3],sep="\t")
  #cat("\n")
}
close(fh)

pattern="\\|\\s+(.\*)\\s+\\["

commapattern="\\|\\s+(.\*)\\,\\s+\\["

getstring=function(text,match,index){
  v=regmatches(text,match)
  vec=v[[1]]
  return(vex[[index]])
}

mymatch=function(pattern,text){
  m=regexec(patter,text,perl=T)
  return(m)
}

ifmatched=function(m){
  if(m[[1]][1]==-1){
    return(FALSE)
  }else{
    return(TRUE)
  }
}

v=rnorm(100)
median(v)
mean(v)
hist(v)

R画图

用最简单的图形命令来画图,在画图之前还穿插了Git的介绍(我会在后续文章中介绍)。


data("airquality")
mean(airquality$Ozone)
mean(airquality$Ozone[is.na(airquality$Ozone)])

d=airquality$Ozone
index=is.na(airquality$Ozone)
clean_idx=!index
clean_value=d[clean_index]
mean(clean_value)
mean(airquality$Ozone,na.rm=T)
summary(airquality)


plot(airquality$Temp,airquality$Ozone, xlab="Temp", ylab="", main="Temp vs Ozone")

plot(airquality$Temp,airquality$Ozone,type="n",xlab="",ylab="",axes=F)
points(airquality$Temp,airquality$Ozone,pch=16)
axis(1)
axis(2)
box()
title(main="Temp vs Ozone",xlab="Temp",ylab="Ozone")
legend("topleft",legend="Some points",pch=16)

par(mfrow=c(1,2))
plot(airquality$Temp,airquality$Ozone, xlab="Temp", ylab="", main="Temp vs Ozone")
plot(airquality$Temp,airquality$Ozone, xlab="Temp", ylab="", main="Temp vs Ozone")


课程简介

课时8小时(2天),介绍linux。

第一天

首先呢,第一节课就讲了怎么用虚拟机载入已经制作好的虚拟镜像。 第二节是历史课,先对自己做了自我介绍,做实验的,博士后开始转行搞生物信息,纯自学成才(还学的这么厉害,我非常佩服)。 先从什么是window/unix/linux谈起,介绍了这几个的区别和历史进程,介绍了linux的读音,mac的OS和linux的有些类似。 然后开始login用户,学习简单的命令。 第三节课用到了挂在的本地目录地址,在目录下写了第一个shell脚本(helloworld.sh)。 由于安装的是极简版的CentOS,所以先从下载firefox和gedit以及gedit插件开始,然后还修改了gedit的首选项。 shebang line别忘了写。 可能还顺便讲了一下C语言的helloworld写法。

第二天

开始讲更多的linux命令,和ncbi数据库。中间还扯到了Joseph Lister Hill捐助了好多钱用来建立数据库(是哪个我忘了)。

虚拟机载入

首先,这个课程已经提前为大家制作了虚拟镜像,只要下载VirtualBox。

第一步,新建一个虚拟机。根据已有镜像文件CentOS-64bit选择linux和RedHat64bit,给虚拟机起个名字,一路next,但最后要选择加载已有镜像。

img

img

第二步,在镜像创建之后,选择设置,在里面将虚拟机和外部电脑桌面的互动打开。

img

第三步,选择本地磁盘存储,一定要勾选自动挂载(automount),这样在/media/文件夹下可以出现你的本地磁盘目录。例如:/media/sf_CB2/。 在folder path这一项中,选择other,就可以找到本地文件夹了。

img

设置完之后可以运行start,开始玩虚拟机!

关机时可以选择保存现有的操作,这样就不用重复启动系统,直接点开始,就能进入桌面。

img

Linux 命令

CB2讲述的linux相关内容有:

  • –help 查看帮助 ls--help
  • man 查看详细帮助 man ls
  • info 查看详细信息 info ls
  • apropos 在一些特定的包含系统命令的简短描述的数据库文件里查找关键字,然后把结果送到标准输出 apropos ls
  • 重定向和管道 重定向>写入文件;>>追加写入文件; cat abc.txt | less 管道符号|用于链接两个命令,前一个命令的输出将作为后一个命令的输入。
  • cat 标准输出(打印到屏幕)文件内容 cat abc.txt
  • more 类似于cat,more以分页的显示形式显示内容 more abc.txt
  • wc 统计文档的行数,单词数,字节数 wc abc.txt
  • head 用来显示文档的开头内容,标准输出(打印到屏幕) head abc.txt
  • tail 用来显示文档的结尾内容,标准输出(打印到屏幕) tail abc.txt
  • alias
  • 路径(pwd, ., ..) pwd可以显示当前路径
  • 环境变量 (export)
  • 执行程序
  • echo 终端打印输出 echo "hello world"
  • 打包和解压缩(tar, gz, bzip2, xz -_,- 一下子讲这么多,记得过来就怪了
  • cut 剪切文件的某一列或者某几列 cut -f 1 abc.txt (剪切第一列,并输入第一列到屏幕)
  • sort 文件按某列排序 cat abc.txt | sort
  • uniq 文件按某列去除冗余 cat abc.txt | sort | uniq
  • wget 下载网上的文件
  • grep 查找关键词 grep "hello" abc.txt
  • tr 对来自标准输入的字符进行替换、压缩和删除。
  • find 查找文件
  • rsync 海量文件的快速准确备份,复制,删除
  • unison 文件同步工具
  • ls 目录文件列表
  • cp 拷贝
  • rm 删除
  • mv 移动
  • cd 进入目录 cd /home/abc
  • mkdir 创建文件夹
  • yum/dnf/apt-get 下载软件包的命令(ubuntu用apt-get)
  • chmod 改变文件的权限

Hello World! 脚本


#!/bin/bash

echo "Hello World!"

课程用到的命令

所有下面的工作最好都在挂载的外部文件夹中操作

cd /media/sf_CB2

下载软件

yum install firefox
yum downgrade http://vault.centos.org/7.1.1503/os/x86_64/Packages/pygobject3-3.8.2-6.el7.x86_64.rpm http://vault.centos.org/7.1.1503/os/x86_64/Packages/pygobject3-base-3.8.2-6.el7.x86_64.rpm
yum install gedit gedit-plugins

下载蛋白质氨基酸序列

wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/proteomes/9606.tsv.gz

查找上述氨基酸序列文件中第七列名称中重复次数最高的名称

zcat 9606.tsv.gz | tail -n+4 | cut -f 7 | sort | uniq -c | sort -gr

下载ncib斑马鱼数据

wget -r -A.faa ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Drosophila_melanogaster/RELEASE_5_48/

将斑马鱼数据复制到当前目录下的faa文件夹中

find ftp.ncbi.nlm.nih.gov/ -name "*.faa" -exec cp {} ./faa \;

安装tree,查看目录结构

yum install tree

tree ./

将所有.faa文件合并成一个dm.faa文件

cat ./faa/*.faa > dm.faa 

计算dm.faa文件中的氨基酸序列条数

cat dm.faa | grep ">" | wc -l

计算dm.faa文件中的氨基酸数量

cat dm.faa | grep -v ">" | tr -d "\n" | wc -c

作业

Problem 1

wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/proteomes/9606.tsv.gz
zcat 9606.tsv.gz | tail -n+4 | cut -f 7 | sort | uniq -c | sort -gr

Problem 2

wget -r -A.faa ftp://ftp.ncbi.nih.gov/genomes/archive/old_refseq/Bacteria/Yersinia_pestis*

Problem 3

find ftp.ncbi.nih.gov/genomes/archive/old_refseq/Bacteria/ -name "*.faa" -exec cp {} ./ypest \;

cat ./ypest/*.faa >  ypest.faa

cat ypest.faa| grep ">" | wc -l

Problem 4

下载

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.faa

写脚本计算

提示:

Performing Math Calculation in Bash

shell float number in expr

How to use expr on float?

#!/bin/bash

faa=$1

c=`cat $faa | grep ">" | wc -l`

l=`cat $faa| grep -v ">" | tr -d "\n" | wc -c`

result1=`expr $l/$c`
result2=$(expr "$l"/"$c")
result3=$(($l/$c))
result4=$(awk 'BEGIN{print $l/$c}')
result5=$(echo "scale=6;$l/$c"|bc)

echo $result1
echo $result2
echo $result3
echo $result4
echo $result5


好长时间没有写课程相关的东西了,学过好几门生物信息学课程,从中可以学到不同老师的教学风格。现在回首,看看研究生时期的生物信息学II(动物所韩老师 -_-b 第一次错写成春雨了,自己都吓了一跳),仅一次内容,介绍了LAMP的初级核心命令,照顾了不同平台尤其是windows平台的同学(在课上使用Cygwin和PuTTY)。

这学期学Bioinformatics computing,每次课程内容很多,连上4个小时,让非计算背景的人无从下手(尤其是连windows安装软件都困难的360用户们对,我就是在喷国内的流氓软件)。 虽然内容极为丰富使用,但没有照顾到非linux用户,前两次课程效果不好。

先对比一下之前的生物信息学II,讲述了linux的相关内容有:

  • PuTTY
  • passwd
  • exit
  • pwd
  • ls
  • man
  • cp
  • vi
  • more
  • cd
  • mkdir
  • su
  • chmod
  • chgrp
  • usermod
  • cat
  • rm
  • find
  • 执行程序
  • 如何进入www目录
  • mysql
  • php (简单带过)
  • alias

CB2讲述的linux相关内容有:

  • –help
  • man
  • info
  • apropos
  • 重定向和管道
  • cat
  • more
  • wc
  • head
  • tail
  • alias
  • 路径(pwd, ., ..)
  • 环境变量 (export)
  • 执行程序
  • echo
  • 打包和解压缩(tar, gz, bzip2, xz -_,- 一下子讲这么多,记得过来就怪了
  • cut
  • sort
  • uniq
  • wget
  • grep
  • tr
  • find
  • rsync
  • unison
  • ls
  • cp
  • rm
  • mv
  • cd
  • mkdir

两个列表凑一凑,就是最全的常用命令集合了。

如果让我讲这个课程的话,如果非要统一使用linux,那就绝对不使用虚拟机,直接用u盘安装可写入型的linux系统安装盘。每次上课,插上u盘换系统,直接使用,可以更新安装软件,所有文件都保存在挂载的电脑磁盘上。可以省掉不少麻烦。

这次还来回换了2次系统,去年上课用fedora,今年上课,一开始让我装32位xubuntu,结果上课当天换成了64位的CentOS,结果我帮助的人在建立虚拟机时我都给人家装错了(RedHat/Ubuntu)_-_

还有最开始的学的生物信息学I,这个怎么评价呢?院士的课怎么学都是好(听老师侃侃而谈很有意思),算法部分还是很有用的,虽然都是课下自学,当时在礼堂课上根本就听不懂,我问的问题竟然由于声音太小,在视频里什么都听不见orz。不知道他老人家现在还开不开课了。

不多说了,赶紧写总结,希望这次BioComp的课程学生对老师的评价不要太低。-,-