2020-01-09 使用Unix进行数据分析的演示

操作酵母菌基因组数据库（SGD）中的SGD_features.tab文件，文件描述：“染色体特征，坐标，基因描述。”

示例数据来源链接
https://downloads.yeastgenome.org/curation/chromosomal_feature/SGD_features.tab
http://data.biostarhandbook.com/data/SGD_features.tab
http://data.biostarhandbook.com/data/SGD_features.README
wget或curl从指定链接下载

# Create a work directory for a unix tutorial work.
mkdir unixtut
# Switch to the new directory you will work in.
cd unixtut
# Get the data from SGD.
wget http://data.biostarhandbook.com/data/SGD_features.tab
# Also get the README file.
wget http://data.biostarhandbook.com/data/SGD_features.README
# Quick check. Where are you and what have you got?
pwd
ls

来源于Unix的命令可以使用man打开帮助手册查看可用的命令行选项
其他软件来源的命令可以尝试-h -help --help等获取帮助

• 命令行选项格式（后称作flags formats）

1. -和字母组合，如：-h
2. --和单词组合，如：--help
3. 某些生信工具会将两者自由组合，如：-genome
4. 有些flag后会有额外的参数，如：-p和-p 100

cat 看README文件

less SGD_features.README

cat 通过管道 | 传递命令给 wc -l， 显示数据有多少行

cat SGD_features.tab | wc -l 
#如下可以执行同样的操作
wc -l SGD_features.tab 
#但更推荐首先使用cat打开数据流，这样设计、构建和理解复杂的管道会更容易。

cat 查看tab文件

cat SGD_features.tab #大量内容从眼前闪过，根本看不清
#head 显示数据的头部
cat SGD_features.tab | head

grep 匹配数据中包含基因YAL060W相关的行

cat SGD_features.tab | grep YAL060W
#显示匹配到的行，高亮匹配词，显示头部
cat SGD_features.tab | grep YAL060W --color=always | head
#这里 --color=always在pattern “YAL060W”有点奇怪，改变顺序也不影响输出结果

grep -v 显示匹配不上的行数

cat SGD_features.tab | grep -v Dubious | wc -l

> 将匹配到的内容写入新的文件

cat SGD_features.tab | grep YAL060W > match.tab
ls 
#这时候应该多了一个match.tab文件

cut -f 取出某一列

cat SGD_features.tab | cut -f 2 | head #第二列为feature type
cat SGD_features.tab | cut -f 2 | grep ORF | head #匹配出feature type为ORF的行，ORF代表可转录的基因
#gene的个数
cat SGD_features.tab | cut -f 2 | grep ORF | wc -l

cut -f 可以一次选取多个不连续的列，并进行操作

cat SGD_features.tab | cut -f 2,3,4 | grep ORF | head 
cat SGD_features.tab | cut -f 2,3,4 | grep ORF | grep -v Dubious | wc -l

探索该数据集中的feature type

#备份出1个新文件，在新文件中操作
cat SGD_features.tab | cut -f 2 > types.txt
#sort 对行进行排序
cat types.txt | sort | head
#uniq 去除重复的行
cat types.txt | sort | uniq | head
#uniq -c 同时给出重复的行的数量
cat types.txt | sort | uniq -c | head
#sort -n 按数字排序
cat types.txt | sort | uniq -c | sort -n | head

#sort -r 反向排序
cat types.txt | sort | uniq -c | sort -rn | head

#发现CDS（蛋白质编码序列）最多7074个,
#其次为ORF（开放阅读框）6604个,
#后面依次为非编码外显子、长末端重复序列、内含子、ARS（自主复制序列）...

sort | uniq -c | sort -rn 是一个非常有用的命令组合

因为uniq只对临近的重复的行起作用，因此要sort对行预先进行排序
有很多有用的命令组合，因此要多学习别人的代码，发现那些你需要的。