2.2.Practice Guide

本章内容主要是学会用Linux的一些简单编程方法去查看GTF/GFF基因组注释文件的基本信息,并学会对文件中数据进行提取,利用提取到的数据计算特定feature(例如计算基因积累长度等)。

awk & sed

在linux中,一般会自带一些awk、sed等类似于perl和python之类的简单语言,可以直接在命令行中使用,执行一些简单的(例如替换字符)操作。

第一次接触linux的同学可以忽略或推迟下面教程中有“awk”或者“sed”出现的命令行的学习。

I. 准备:了解基因组注释文件

顺利掌握以下Linux操作的前提是,必须熟悉GTF/GFF文件的注释信息,特别是要每一列对应的内容是什么。

1. 获取文件

获取方法

2. gff/gtf文件格式和tab分隔符

  • GFF全称为general feature format,这种格式主要是用来注释基因组, 9列数据组成。

Example:

chr22 TeleGene enhancer 10000000 10001000 500 + . touch1
chr22 TeleGene promoter 10010000 10010100 900 + . touch1
chr22 TeleGene promoter 10020000 10025000 800 - . touch2
  • GTF全称为gene transfer format,主要是用来对基因进行注释,前8列和gff相同的,第9列有一些小的差别。

以下为每一列的对应信息:

Tips: 为了便于处理,数据文件建议都做成矩阵(matrix)的形式,也就是行列清晰。列和列之间一般建议用tab分隔符(键盘上的tab按键)分开,而不是一个空格键分开。

更多解释:http://www.genome.ucsc.edu/FAQ/FAQformat.html

II. 开始正式练习:Linux命令练习

操作要点:

请在操作前和操作后分别阅读一下,仔细体会。

  • 要点1. Linux命令行格式通常写法如下:

命令(空格)选项(空格)参数1 参数2...

mv -f folder1 folder2 #实际使用时将folder1替换为需要移动的文件夹,folder2替换为希望移动到的位置。

注意:命令、选项、参数之间一定要用空格来区分!

  • 要点2. 两种表达方式:短格式 vs 长格式。

    • ① 短格式的命令选项:用一个 - 和一个单个英文字母表示, 如 -a

    • ② 长格式的命令选项:用两个 - 和一个英文单词表示, 如 --help

ls -hls --help 或者 ls -als --all 所起的作用都是相同的。

  • 要点3. cd——进入工作目录

cd #cd后面为空时,进入默认家目录
cd ~/linux #工作目录名称,这里为本章工作目录 ~/linux,TAB键可进行名称自动补全,推荐经常使用

一般的程序后面都要输入文件位置和名称,告知程序输入和输出是什么:

./filename 指当前目录下的文件 ../filename 指上一级目录下的文件 ../../filename 指上两级目录下的文件.

  • 要点4. "|"是管道命令操作符,它可以将左边命令传出的正确输出信息(standard output)作为右边命令的标准输入(standard input)。

  • 要点5. 建议在对*.gtf文件执行的一些命令行Inputs末尾加上| head -n或者| tail -n,然后Outputs会自动显示文件前n行或者后n行;否则,屏幕会被刷屏。

  • 要点6. 星字符:*可以代表任何字符,称之为wildcard。

重点和难点:

awkcatcutgrepwc的用法其参数的用法是本章学习的,也是主要的homework之一。

step0.准备: 解压缩.gtf的文件

cd ~/linux
gunzip 1.gtf.gz
ls # check if 1.gtf.gz has been unzipped to 1.gtf

step1.查看文件基本信息

尝试输入以下命令,分别查看1.gtf文件的开头、结尾、文件的大小、行数等基本信息。

cat 1.gtf | head #显示1.gtf文件前10行
cat 1.gtf | tail #显示1.gtf文件后10行
cat 1.gtf | head -15 #显示1.gtf文件前15行(输入值15可以用其他整数替代)
ls -lh 1.gtf #显示1.gtf文件的大小
wc -l 1.gtf #统计1.gtf文件行数
grep -v "#" 1.gtf | grep -v '^$' | wc -l #用grep -v排除commend line(以#开头的部分)以及无意义的空白行
awk '{if($0!=" ") print}' 1.gtf | head -10 #过滤空行,显示前十行结果。
grep -v '^#' 1.gtf |head -5 #过滤开头注释行commend line(以#开头的部分)并显示前5行

step2.数据提取

首次尝试,先复制以下命令,分别提取1.gtf文件的特定列、行等数据信息;观察输出结果,然后建议尝试修改以下命令中的参数,进行更多的练习。

2.1 筛选特定的列

#选取1-3列的数据(以下两种命令都可以)
cat 1.gtf | awk ' { print $1, $2, $3 } ' | head
cat 1.gtf | cut -f 1,2,3 | head
#Eg.例如我只需要GTF文件的第1,34,5列也就是chr,feature,start,end。
cut -f 1,3,4,5 1.gtf | head

2.2 筛选特定的行

# 假设我们想要提取第三列是gene的行,并且只显示第1,3,9这几列信息。
cat 1.gtf | awk '$3 =="gene" { print $1, $3, $9 } ' |head

Step3.提取和计算特定的feature

这一阶段是在学会step2的基础上,进一步的学习。首次尝试,先复制以下命令,观察输出结果,然后建议尝试修改以下命令中的参数,进行更多的练习。

3.1 提取并统计featrue类型

grep -v '^#' 1.gtf |awk '{print $3}'| sort | uniq -c #提取并计数有多少类feature

3.2 计算特定feature特征长度

#*第5列的数值减去第4列的数值后+1,即得到特征feature的长度
cat 1.gtf | awk ' { print $3, $5-$4 + 1 } ' | head
# 计算所有gene的累积长度
cat 1.gtf | awk '$3 =="gene" { len=$5-$4 + 1; size += len; print "Size:", size } ' | tail -n 1
#计算所有CDS的累积长度
cat 1.gtf | awk '$3 =="CDS" { len=$5-$4 + 1; size += len; print "Size:", size } ' | tail -n 1
#计算1号染色体cds的平均长度
awk 'BEGIN {s = 0;line = 0 } ;$3 =="CDS" && $1 =="I" { s += ($5 - $4);line += 1}; END {print "mean=" s/line}' 1.gtf

3.3 分离并提取基因名字

#从gtf文件中分离提取基因名字,并计算其长度
cat 1.gtf | awk '$3 == "gene"{split($10,x,";");name = x[1];gsub("\"", "", name);print name,$5-$4+1}' | head

step4.提取数据并存入新文件

这一阶段主要是学会提取数据并存入新文件,例如,寻找长度最长的3个exon, 汇报其长度

这里介绍两种方法。

第一种是直接提取并计算最长 3 个 exon, 汇报其长度,存入 .txt 文件;

第二种方法是写一个可执行文件 run.sh,寻找长度最长的 3 个 exon,汇报其长度。

4.1 提取数据存入txt文件示范

输入命令如下,则可将结果存入新文件 1.txt,这里提取并存入的命令是>

grep exon 1.gtf | awk '{print $5-$4+1}' | sort -n | tail -3 > 1.txt

然后输入命令 less 1.txt 或者 vi 1.txt 则可进入vi一般模式界面显示输出结果,例如:

vi简单使用教程详见Tips,此时,在英文输入法状态下按:q:wq可以退回到终端shell窗口。

在输入less查看文件时,也可以使用q退出查看模式。

1.txt 拷贝到/home/test/share,就可以在本地 (~/Desktop/bioinfo_tsinghua_share)查看 1.txt 文件。

4.2 可执行文件编辑示范

第一步,输入命令,进入vi编辑界面。

vi run.sh

第二步,按 i 键切换至insert模式后,写下rush.sh的文件内容如下:

#!/bin/bash
grep exon *.gtf | awk '{print $5-$4+1}' | sort -n | tail -3

(第一行语句一般用来声明这个脚本使用的shell名称,# 后的语句可作为批注,在执行时会被忽略)

第三步,按 Escctrl+[ 切回普通模式,输入 :wq 退出vi编辑器,在命令行后键入:

chmod +x run.sh
./run.sh

输出如下所示,与 1.txt 的内容一致。

III. Homework

  1. 解释1.gtf文件中第4、5列代表什么,exon长度应该是$5-$4+1还是$5-$4?

  2. 列出1.gtf文件中 XI 号染色体上的后 10 个 CDS (按照每个CDS终止位置的基因组坐标进行sort)。

  3. 统计 IV 号染色体上各类 feature (1.gtf文件的第3列,有些注释文件中还应同时考虑第2列) 的数目,并按升序排列。

  • 作业格式:提交word/md/txt/sh文件均可

  • 作业解释:第2,3题要求给出结果,也可以附上使用的命令

  • 注: 第一次接触linux的同学可以忽略或推迟上述教程中有“awk”或者“sed”出现的命令行的学习; 只有一个很简单的语句最好掌握一下:awk '{print $5-$4+1}' #输出第5行的数减去第4列的数后加1.

本章参考文献:

本章主要参考以下几篇生信数据文章: ​https://www.jianshu.com/p/48b5a0972301​ ​https://blog.csdn.net/sinat_38163598/article/details/72851239​ ​https://zhuanlan.zhihu.com/p/36065699​ ​https://gist.github.com/sp00nman/10372555​ ​https://www.jianshu.com/p/7af624409dcd