慢慢看毡庆，憋著急坑赡！很有用！

前言：

首先呢么抗，在你的Linux系統(tǒng)中新建一個(gè)文件毅否，Thanos.txt(紫薯俠賜予你力量)

@敲入這一段，待會(huì)兒肯定給你解釋清楚啥意思
echo -e "Hello world\nI love bioinformatics\n\nWelcome to bioinfoplanet\nAnd I will be happy here\n\nWelcome come to our Wechat public account\nbioinfoplanet" > Thanos.txt

然后你會(huì)發(fā)現(xiàn)多了一個(gè)紫薯俠文件蝇刀，cat Thanos.txt 看是不是輸出了幾行字螟加？真棒??你已經(jīng)開始了第一步

解釋一下，echo是打印到終端，你可以理解為打印機(jī)捆探，我們賦予它內(nèi)容然爆，他就給我們打印出來(lái)，到Thanos.txt這個(gè)本子上黍图。

• -e ： 內(nèi)容中有特殊字符就需要加（比如上面??的換行符\n）
• " "中就是要打印的內(nèi)容
• > 輸出到哪里

利用這個(gè)文件曾雕，你可以繼續(xù)向下學(xué)習(xí)。當(dāng)然如果你認(rèn)為內(nèi)容太少助被，你可以vi Thanos.txt 進(jìn)入vi編輯器輸入i編輯

文本查看：

查看命令

1. head, tail 查看文件頭尾

   -n 查看指定行 (默認(rèn)值為10行)

2. cat 將文件全部打出來(lái)

3. less 逐頁(yè)顯示文本

   -S 規(guī)則輸出

   -N顯示行號(hào)

4. less界面中可以移動(dòng)光標(biāo)剖张，搜索關(guān)鍵字

less可以進(jìn)行debug，查看中間輸出結(jié)果
例如有一個(gè)管道程序sh1 test.txt | sh2 | sh3 > output.txt
【其中sh1-3都是命令】
在運(yùn)行之前揩环，我們可以看看每一步是否正確搔弄，這樣使用：
sh1 test.txt | less
sh1 test.txt | sh2 | less
sh1 test.txt | sh2 | sh3 | less

1. 查看行數(shù)：wc -l 輸出行數(shù) （但是空行也會(huì)被算入）=> 解決： grep -c "[^ \\n\\t]" file
2. 去除注釋信息（元數(shù)據(jù)）后查看行數(shù)
   一般基因組注釋GTF文件開始都有幾行的注釋文件，需要去除以后再查看
   grep -v '^@' *.gtf | wc -l
   這樣做的好處就是接下來(lái)還能用awk直接處理列或者計(jì)算列數(shù)
   grep -v '^@' *.gtf | awk -F "\t" '{print NF};exit'
3. 查看列數(shù)【以tab分割為例】：
   awk -F "\t" '{print NF ; exit}' *.bed

編輯器命令：

q：退出

g：第一行 G：最后一行

/<pattern> 向下??搜索 丰滑？<pattern> 向上??搜索

n向后匹配 N向前匹配

文件輸出：

echo -n 不換行輸出 -e 處理特殊字符

echo -e "nihao\nworld"                  echo -n "nihao\nworld"                      
=> nihao                                => nihao\nworld
=> world

>覆蓋原文件 >>添加到原文件底部

正則表達(dá)式regular expression：

它的語(yǔ)法結(jié)構(gòu)有兩套系統(tǒng)組成顾犹，元字符（metacharacters） + 普通字符

那么問(wèn)題來(lái)了：什么是元字符呢？
A： 元字符是表達(dá)式的結(jié)構(gòu)吨枉，骨架蹦渣；比如，“我愛(ài)生信星球”這個(gè)句子中的主謂賓都是固定的貌亭，也就是它的元字符是固定的；普通字符可以不同认臊，也就是會(huì)有各種語(yǔ)言版本的這句話圃庭，來(lái)表達(dá)我對(duì)生信星球的愛(ài)

元字符主要由以下字符組成：^ $ . [] {} - ? + () | \

()與[]的差別：()為多選儒旬，[]為單選, 對(duì)比下就知道

cat SRR519926_1.fastq | egrep 'TA(A+)TA' --color=always | head
@ 匹配了多種模式

re_1.png

cat SRR519926_1.fastq | egrep 'TA[A+]TA' --color=always | head
@ 僅匹配了一種模式

re_2.png

(e)grep：號(hào)稱“Almost～最快文本搜索”

意思是 global search regular expression(RE) and print out the line
egrep 是grep的拓展模式，支持的元字符較多

生而為搜——它會(huì)在每一行搜索匹配內(nèi)容

-c 顯示有多少行被匹配到（count）

-v 過(guò)濾掉某些格式的行

-w 完全匹配【比如想刪除匹配到abc的帖族，不使用-w可能abc1和abcd都要被刪除】

-o 只打印匹配到的內(nèi)容

-i 忽略匹配字段和匹配內(nèi)容的大小寫

-A/B n： 輸出匹配內(nèi)容前/后 n行

--color=always/auto: 始終/自動(dòng)高亮顯示搜索字段

• 搜索特定信息

  @ 搜索特定基因信息(在擬南芥gff文件中查找栈源，下載地址見下方??)
  grep "AT1G01680" *.gff | head -n5
  

grep_1.png

• 排除特定信息

  @ 還是搜索這個(gè)基因，但排除其中feature項(xiàng)的protein
  grep "AT1G01680" *.gff | grep -v "protein" | head -n5
  

grep_2.png

• 查找某段序列并輸出上下文

  @ grep -A n顯示后面n行竖般； -B n顯示前面n行
  @ 下載基因組文件wget http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas
  
  grep -A 2 "CAAATTGAATTAAG" *.fas
  => 得到下面??的結(jié)果
  => -B的用法你可以試一下
  

grep_3.png

• 查找特定序列并計(jì)算出現(xiàn)了幾次

  grep -c "CAAATTGAATTAAG" *.fas 
  或者： grep "CAAATTGAATTAAG" *.fas | wc -l
  @ 如果單純輸出 就用grep -o
  

• 精準(zhǔn)匹配某個(gè)基因

  @ 搭配正則表達(dá)式
  @ 比如要匹配開頭為AT1G250甚垦，結(jié)尾為3的基因名
  grep "AT1G250.*3$" *.gff
  => 得到了兩條AT1G25053 & AT1G25083
  

grep_4.png

• 并非完美～

  雖然grep很強(qiáng)大，但是他也并非十全十美，他的一個(gè)不足之處就在于艰亮，對(duì)于存在換行符的字符串，它會(huì)搜不到迄埃。例如我們想找TAIR10_chr_all.fas 中95-100之內(nèi)的‘CCACT’ 堿基

  先看一下tail -n100 TAIR10_chr_all.fas | head -n5

grep_5.png

用肉眼就能看到'CCACT'疗韵，但是如果用grep搜索 是沒(méi)有結(jié)果的
tail -n100 TAIR10_chr_all.fas | head -n5 | egrep 'CCACT'

其實(shí)稍加改變伶棒，用grep還是能實(shí)現(xiàn)，就是不那么優(yōu)雅：
tail -n100 TAIR10_chr_all.fas | head -n5 | tr -d '\n' | egrep 'CCACT'
tr -d（delete） 是刪除特定字段

所以眯搭，有沒(méi)有什么更快捷的辦法呢窥翩？
有的～可以使用emboss套件下的dreg, 它是針對(duì)核酸鳞仙；如果是氨基酸序列和翻譯后的蛋白序列寇蚊，使用preg

tail -n 1000 chr22.fa | head -n 5 | dreg -filter -pattern 'TAATA'

sed:

流編輯器～就像水流一樣，按行從數(shù)據(jù)讀取--執(zhí)行命令--顯示結(jié)果棍好，一氣呵成
整個(gè)過(guò)程在一個(gè)位于內(nèi)存中叫Pattern Space模式空間中進(jìn)行仗岸，因此不會(huì)改變?cè)嘉募?nèi)容

使用方法：
sed [options] 'Adress Command1;Command2;...' file...

[options]包括：
-n: 只顯示經(jīng)過(guò)sed處理的行，保存在模式空間的未處理行不顯示【常用】
-i：直接修改讀取內(nèi)容借笙，不輸出【慎用扒怖，他會(huì)修改源文件】
-r：使用拓展正則表達(dá)式

Adress工作范圍：
可以指定為按行號(hào)工作或按過(guò)濾條件工作

• 按行號(hào)：使用單獨(dú)的數(shù)字表示某一行；兩個(gè)數(shù)字并逗號(hào)分隔表示范圍业稼；
  n,m+ 表示從n向下m行盗痒；n~m表示從n開始的每m行

  比如我要輸出偶數(shù)行: sed -n “2~2 p" Thanos.txt

• 按匹配條件：'/pattern/command' 在包括pattern的行中執(zhí)行命令

  比如要打印包含bio的行：sed -n '/bio/p' Thanos.txt
  比如打印包含bio與包含planet之間的所有行：
  sed -n '/bio/, /planet/p' Thanos.txt
  還可以與行號(hào)連用：sed -n '/bio/, +6 p'Thanos.txt 打印包含bio行以及下面6行

Command命令[多個(gè)命令用分號(hào)分隔]：
p：復(fù)制模式空間的內(nèi)容，一般與-n連用【否則會(huì)一次輸出兩次】
d：刪除【按行號(hào)或匹配條件】

刪除第5行后面所有的行
sed '5,$d'

i：插入

例如低散，要在文件開頭插入一行俯邓，三列name、length谦纱、feature
sed '1i name\tlength\tfeature' Thanos.txt

a：追加

例如看成，要在末尾追加內(nèi)容【表示最后一行】 `sed 'a bioinfoplanet' Thanos.txt`

c：替換【與a用法相似】
n：匹配的行向下移一行進(jìn)行操作

例如：seq 6 | sed '2{n;d}' 輸出的就是除了3以外的數(shù)字
再如：seq 6 |sed 'n;d' 輸出奇數(shù)
seq 6 | sed -n 'n;p' 輸出偶數(shù)

!：反向執(zhí)行

sed '/bio/!d' Thanos.txt 效果等于 sed -n '/bio/p' Thanos.txt

=：打印行號(hào)

sed '/bio/!d;=' Thanos.txt

s：替換【支持正則】
s/pattern/replace/flags

其中pattern是支持正則的
flags包括：n（替換第n個(gè)匹配）；g（全局匹配）跨嘉；p（與-n搭配川慌，輸出修改的行）吃嘿；i（忽略大小寫）；w（保存修改的行到一個(gè)新文件）

一種特殊情況梦重，還比較常用兑燥，就是如果替換/怎么辦？
方法一：使用轉(zhuǎn)義符\/
方法二 ：用@ | ! ^ 替換
例如：要更改當(dāng)前目錄中的部分內(nèi)容
pwd | sed 's@/home/tmp/bio@/home/tmp/bioinfo@'

還有一種特殊情況琴拧，就是有時(shí)候不想替換掉降瞳，只是想把pattern的這部分內(nèi)容與replace的內(nèi)容一起輸出來(lái)：
sed 's/[[:upper:]]/ word = &/' Thanos.txt
意思就是：將Thanos.txt中的大寫字母放到&的位置，輸出格式就是：word = 大寫字母蚓胸≌跫ⅲ【特殊字符&用來(lái)存儲(chǔ)pattern中的內(nèi)容】

實(shí)用的sed單行命令：

@ 刪除空行
sed '/^$/d' 

@ 每一行下增加一行空行
sed G

@ 每三行增加一行空白行
sed '0~3G'

@ 在匹配的pattern后面加入一行空白行
sed '/pattern/G'

@ cat的功能實(shí)現(xiàn)
sed ''

@ head的功能實(shí)現(xiàn)
sed '3 q' =》輸出前三行

@ tee功能實(shí)現(xiàn)
sed -n 'p; w newfile'

@ grep功能實(shí)現(xiàn)
sed -n '/pattern/p'

@ grep -v功能實(shí)現(xiàn)
sed -n '/pattern/p!'

@ 計(jì)算行數(shù)
sed -n '$='

@ 多個(gè)內(nèi)容同時(shí)替換(例如將1、2沛膳、3替換成4)
sed 's/1\|2\|3/4/'

@ 顯示包含“haha”扔枫、“xixi”、“yeah”的行
sed '/haha/d!; /xixi/d!; /yeah/d!'

awk：號(hào)稱“Almost～最強(qiáng)文本操作”

“awk w(o/a)rd”—— 讓看似繁復(fù)無(wú)章的文字遇見它就秒變尷尬??

工作職責(zé)：主打行【如果一個(gè)文本文檔是一張表格锹安，每一行代表一個(gè)記錄短荐，每一列代表域，awk就是處理記錄專用】

工作流程： 先逐行讀取并記錄叹哭，將行信息整合入$0 => 指定分隔符（默認(rèn)空格）忍宋，分割成各個(gè)列$1,$2,$3...；

再執(zhí)行awk 'pattern1 {command1}; pattern2 {command2}...'

??：情況一：如果缺失pattern【也就是利用下面幾種運(yùn)算進(jìn)行的匹配】风罩，直接執(zhí)行{command}糠排；
情況二：缺失{command}【在{}中輸入的命令，如print】超升，執(zhí)行pattern

現(xiàn)在不太懂也沒(méi)關(guān)系乳讥，通過(guò)實(shí)例練習(xí)你就明白，下面的實(shí)例我會(huì)描述對(duì)應(yīng)哪種情況

• 模擬cat打印行

  awk '{print $0}' Thanos.txt | head -n3 
  @ 【因?yàn)闆](méi)有pattern廓俭，所以直接執(zhí)行打印命令】～屬于情況一
  

• 選擇列打印

  @ 利用.gff文件練習(xí) ～屬于情況一
  awk '{print $1,$3,$4,$5}' *.gff | head -n5
  【當(dāng)然，這里默認(rèn)是空格分割唉工，如果想tab分割呢研乒？】
  awk '{[print $1 "\t" $3 "\t" $4 "\t" $5 "\t"]}' 
  @ 一會(huì)練習(xí)完下面的cut命令，看看有什么發(fā)現(xiàn)淋硝？
  

awk_1.png

• 列重排

  @ 這個(gè)功能相當(dāng)于下面cut的升級(jí)版 ～ 屬于情況一
  @ 可以改變列的順序雹熬，并且可以自定義分隔符
  awk '{print $4","$5":"$1":"$3}' *.gff | head -n5
  @ 隨心所欲，只為你改變谣膳！
  

awk_2.png

• 轉(zhuǎn)換格式

  比如將gff/gtf格式轉(zhuǎn)為bed格式
  【注意兩點(diǎn)：一竿报、轉(zhuǎn)換后的bed格式分隔符為tab；
  二继谚、bed與gff的坐標(biāo)格式不同：0-坐標(biāo)系統(tǒng)有BAM烈菌、BED、BCFv2、PSL芽世；1-坐標(biāo)系統(tǒng)有SAM挚赊、VCF、GFF济瓢、wiggle荠割。需要做出調(diào)整】
  cat test.gtf | awk '{print $1 "\t" $4-1 "\t" $5}' > test.bed

• 如果單有這些簡(jiǎn)單的功能，還算不上Almost最強(qiáng)旺矾，加上下面這些你再試試蔑鹦？

  劉小澤喊你進(jìn)階啦！

  • 邏輯運(yùn)算（<, > , <=, >=, ==, !=）

  • 數(shù)學(xué)運(yùn)算（+箕宙，-嚎朽， *， /扒吁， %）

  • 關(guān)系運(yùn)算（與&&火鼻， 或||， 非５癖馈）

  • 正則運(yùn)算 （實(shí)則就是將你想匹配的東西放在/ /里魁索，然后在它前面加～匹配， !~ 不匹配）

  • • 舉個(gè)小例子1：

      @ 我想匹配在3號(hào)染色體上盼铁，長(zhǎng)度大于1.5k的注釋,看前五行    
      @ ～情況二粗蔚，我們只需要pattern就好
      awk '$5-$4 > 15000 && $1 ~/Chr3/' *.gff | head -n5
      

awk_3.png

  > 然后自己再試試匹配1號(hào)或2號(hào)染色體上，長(zhǎng)度小于1.5k的注釋饶火，默認(rèn)輸出前10行

- 小例子2:

  ```平均值
  @ 我想看一下3號(hào)染色體上編碼區(qū)（CDS）的平均值 ～情況二
  awk 'BEGIN {len=0;line=0}; $1 ~/Chr3/ && $3 ~/CDS/ { len += ( $5 - $4 );line += 1}; END {print "CDS_mean=" len/line}' *.gff
  =>
  CDS_mean=225.019
  【這里 $1 ~/Chr3/ 中的～是為了匹配正則表達(dá)式】
  @ 那么你試試看鹏控，1號(hào)染色體上外顯子exon的平均長(zhǎng)度吧
  @ 此處只是用來(lái)練習(xí)命令，并不表示真實(shí)長(zhǎng)度肤寝。因?yàn)橹丿B的區(qū)域（比如起始位點(diǎn)一樣当辐，終止位點(diǎn)不同）這樣會(huì)被記錄兩次，結(jié)果是不準(zhǔn)確的鲤看；如果要精確統(tǒng)計(jì)缘揪，需要用編程去重以后來(lái)實(shí)現(xiàn)
  ```

• 特殊的變量：

  • （上面說(shuō)的）Num**表示哪一個(gè)字塊（這里再次提醒：**0代表當(dāng)前內(nèi)容）

  • NR代表當(dāng)前所在行號(hào)，想要打印3-4行內(nèi)容：head -4|tail -2义桂，
    就等同于awk 'NR>=3 && NR<=4 {print $0}'`

  • NF表示目前的記錄被分割的字段的數(shù)目找筝，即Number of Field

    例如要輸出文件的列數(shù)【默認(rèn)空格分割，這里設(shè)為tab分割】：
    awk -F “\t” '{print NF; exit}' test.bed

  • FS 指定列分隔符

  • OFS 指定列輸出分隔符

    例如將上游文件的默認(rèn)分隔符\t更換成慷吊， 并輸出3-5列
    cat test.txt | awk '{FS="\t"; OFS=",";}{print $3,$4,$5}'

常用的還有：

cut：提取

• -f 提取指定字段（filed）cut -f 1,3,4,5 GFF3_genes.gff | head
  也可以cut -f 1-5

• -d: 指定分隔符（默認(rèn)\t）

• -c 截取一定范圍的字符袖裕，例如cut -c1-3就是截取三個(gè)字符

  例如，將bed文件第一列的染色體編號(hào)去掉
  awk '{print 1}' *.bed | cut -c4 反之溉瓶，要再添加上的話： `awk '{print1}' *.bed | cut -c4 | awk '{print "chr"$1}'`

uniq: 去重

與sort連用急鳄，-c 在每列旁邊顯示該行重復(fù)出現(xiàn)的次數(shù)
-d ： 只輸出重復(fù)的行

column:格式化輸出

一般cut后的結(jié)果參差不齊谤民，可以用它對(duì)齊輸出結(jié)果，默認(rèn)\t
指定分隔符：以逗號(hào)為例攒岛， -t -s ','
【??：他的使用只是為了好看赖临，不要把它c(diǎn)olumn處理的結(jié)果交給下游繼續(xù)處理，這樣會(huì)讓文本解析速度下降灾锯。??看來(lái)Linux還是更習(xí)慣參差不齊的文本】

sort: 排序

• -k 被排序的部位兢榨，包含start 和stop兩部分, 二者,隔開；每次使用都要加上這兩部分【默認(rèn)排序第一列】
• -n 按數(shù)字大小排序
  sort -k1,1 -k2,2n my.bed 對(duì)第一列首字母按字符排序顺饮，對(duì)第二列按數(shù)值
• -r 降序（默認(rèn)升序）【與-k連用表示對(duì)某一列逆序】
• -t 指定字塊分隔符（默認(rèn)tab）
• -c 檢查是否按照某種方式排過(guò)順序【echo $? 返回0表示執(zhí)行成功】
• -V 排序時(shí)不用ASCII碼方式排列吵聪，就顯示為正常的排序方式
  【比如三個(gè)染色體編號(hào)chr1, chr2, chr11。不加-V排序結(jié)果就是chr1-- chr11 -- chr2; 加了-V就是：chr1 -- chr2 -- chr11】
  【一般數(shù)據(jù)處理過(guò)程中不用這個(gè)操作兼雄，原始的格式處理更快】

join：連接

要求兩個(gè)文件之間必須有共同點(diǎn)吟逝，所以使用join前必須先將文件排序

格式：join -1 <file1_field> -2 <file2_field> <file1> <file2>

共同點(diǎn)通過(guò)-1 、 -2傳遞進(jìn)來(lái)赦肋，比如說(shuō)兩個(gè)文件的某一列有共同點(diǎn)

例如一個(gè)bed文件前三列如下：my.bed
chr1 34 36
chr2 38 40
chr1 25 39
chr3 12 18
一個(gè)染色體長(zhǎng)度文件為：length.txt
chr1 54362
chr2 35613
chr3 46612

join這兩個(gè)文件之前先sort
sort -k1,1 my.bed > sorted.bed
sort -c -k1,1 length.txt
以雙方第一列為共同點(diǎn)拼接：
join -1 1 -2 1 sorted.bed length.txt > with_length.txt
【假如兩個(gè)文件沒(méi)有共同點(diǎn)块攒，比如length.txt少了chr2數(shù)據(jù)，那么join后的文件也缺少chr2數(shù)據(jù)】

-a 選項(xiàng)指定哪個(gè)文件可以不遵循佃乘，在沒(méi)有共同點(diǎn)時(shí)可以單列出來(lái)

舉個(gè)??【很大的那種囱井！】：

首先輸出一下未排序文件，你可以自己下一個(gè)小數(shù)據(jù)(43M)試試, 人類的比較大(1.2G)我就不下了

wget http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff【這是擬南芥的基因注釋文件趣避，當(dāng)練手很好用】

先練習(xí)下之前學(xué)的一些處理：

@0 拿到手?jǐn)?shù)據(jù)庞呕，可以先ll -h *.gff 看一下文件大小
                less -SN *.gff 看一下文本內(nèi)容

@1 去除多余的#注釋行 與 空行， 打印行號(hào)
grep -v "#" *.gff | grep -v "^$" | wc -l
=> 590264行（一般這種使用量大的數(shù)據(jù)是沒(méi)有多余行的程帕，但是人類的有住练。這里只是為了演示用法）

@2 截取gff文件的1-5列，將第二列除去, 輸出默認(rèn)的前10行到test.txt中
cut -f 1,3,4,5 *.gff > test.txt 
cat test.txt | head
=>
Chr1  chromosome  1   30427671
Chr1  gene    3631    5899
Chr1  mRNA    3631    5899
Chr1  protein 3760    5630
Chr1  exon    3631    3913
Chr1  five_prime_UTR  3631    3759
Chr1  CDS 3760    3913
Chr1  exon    3996    4276
Chr1  CDS 3996    4276
Chr1  exon    4486    4605

接下來(lái)對(duì)截取的test.txt進(jìn)行處理

@3 想根據(jù)第二列的feature進(jìn)行排序【注意sort -k的使用愁拭！】
sort -k 2,2 test.txt | head -n5
=>
Chr1  CDS 3760    3913
Chr1  CDS 3996    4276
Chr1  chromosome  1   30427671
Chr1  exon    3631    3913
Chr1  exon    3996    4276

@4 【進(jìn)階】先根據(jù)第一列Chr數(shù)字大小降序排序讲逛，再根據(jù)第二列排序
你是不是試過(guò)了這個(gè)》sort -k 1,1nr -k 2,2 *.txt | head -n5
你會(huì)發(fā)現(xiàn)輸出的結(jié)果第一列還是Chr1開頭，并沒(méi)有降序岭埠，為什么呢妆绞？
=》原因就是，-k 1枫攀，1還是根據(jù)第一個(gè)字段的全部排序，還是根據(jù)Chr1的‘C’進(jìn)行匹配株茶，其實(shí)我們只想用第一個(gè)字段的數(shù)字（也就是第一個(gè)字段的第四個(gè)字符）進(jìn)行匹配
=》如何實(shí)現(xiàn)来涨？其實(shí)sort -k參數(shù) 內(nèi)置了這個(gè)功能。使用-k 1.4启盛，1.4 就是根據(jù)這種特定方式匹配啦
sort -k 1.4,1.4nr -k 2,2 *.txt | head -n5
=> 
Chr5  CDS 10001590    10001736
Chr5  CDS 10004720    10004824
Chr5  CDS 10004720    10004824
Chr5  CDS 10005070    10005255
Chr5  CDS 10005070    10005255
有沒(méi)有很好用蹦掐？技羔！

如果只是想統(tǒng)計(jì)一下整體的feature（第三列）情況，可以這樣：

cut -f 3 *.gff |sort|uniq -c >feature.txt
=>
 197160 CDS
      7 chromosome
 215909 exon
  34621 five_prime_UTR
  28775 gene
  35386 mRNA
   3911 mRNA_TE_gene
    180 miRNA
    480 ncRNA
  35386 protein
    924 pseudogene
   1274 pseudogenic_exon
    926 pseudogenic_transcript
     15 rRNA
     13 snRNA
     71 snoRNA
    689 tRNA
  30634 three_prime_UTR
   3903 transposable_element_gene

想看哪個(gè)feature最多卧抗？沒(méi)問(wèn)題藤滥，一步搞定！

sort -k 1r feature.txt
=>
 215909 exon
 197160 CDS
  35386 protein
  35386 mRNA
  34621 five_prime_UTR
  30634 three_prime_UTR
  28775 gene
   3911 mRNA_TE_gene
   3903 transposable_element_gene
   1274 pseudogenic_exon
    926 pseudogenic_transcript
    924 pseudogene
    689 tRNA
    480 ncRNA
    180 miRNA
     71 snoRNA
     15 rRNA
     13 snRNA
      7 chromosome

【小練習(xí)1：】GENCODE下載人類基因組GRCh38注釋gff3 ,正好練習(xí)下數(shù)據(jù)庫(kù)使用

? 然后統(tǒng)計(jì)人類的基因組 feature

想偷懶社裆？給你個(gè)機(jī)會(huì)～ wget -c ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz

應(yīng)該得到如下結(jié)果：=>

1237914 exon

735618 CDS

203835 transcript

148007 five_prime_UTR

144591 three_prime_UTR

85439 start_codon

77451 stop_codon

58381 gene

119 stop_codon_redefined_as_selenocysteine

【小練習(xí)2:】 相信你能夠完成上面的內(nèi)容拙绊，輸出了結(jié)果。

小前言：我們知道泳秀，Linux基于Unix開發(fā)标沪，肯定要繼承unix的精華，那就是利用小程序的整合去完成大任務(wù)嗜傅。Linux中的管道命令|就是這樣一種體現(xiàn)金句。管道中的數(shù)據(jù)可以不被寫入磁盤，在更高速的內(nèi)存中進(jìn)行處理吕嘀，就像空中綠道违寞，不被道路上擁擠的車流阻攔，你可以在空中自由騎行偶房，肯定效率高很多趁曼，心情也更舒暢～

問(wèn)題來(lái)了：<u>你能否用一行管道命令從讀取文件開始到輸出排序好的feature文件呢？</u>

---

提示：讀取|截取|排序|統(tǒng)計(jì)|統(tǒng)計(jì)后排序