生信分析過程中，會與很多不同格式的文件打交道，除了原始測序數據fastq之外，還需要準備基因組文件fasta格式和基因注釋文件gtf格式。在分析的過程中還會有眾多中間文件的生成，如bed、bed12、sam、bam、wig、bigwig、bedgraph等，生成后我們一般會查看下內容了解文件每一列的含義，以此來決定需要提取哪些有用信息列來進行下一步分析。

插播一個小劇場

老板：“先查看一下bam文件內容。”
小白：嗒嗒嗒敲鍵盤。

$ less ehbio.bam "ehbio.bam" may be a binaryfile. ?See it anyway?

小白：“哎呀，錯了，應該這樣。” ?嗒嗒嗒敲鍵盤。

$ samtools view ehbio.bam # 回車一敲，災難，電腦要卡死，趕緊按`control ?+ c` $ samtools view ehbio.bam | less ?# 這下終于可以查看了

老板：“你逗我呢……”（不失禮貌的批評）

剛接觸生信分析的小白們這種尷尬的事情時有發生，為了幫助大家梳理這些剪不斷理還亂的文件，本文以分析流程為主線，介紹各文件的格式以及有哪些常用命令來查看或處理它們。

1. 測序數據FASTQ文件

1）文件用途：樣品測序返回的數據一般存儲為fastq文件，通常是壓縮文件filename.fq.gz的格式，節省存儲空間和傳輸時間。NGS基礎 - FASTQ格式解釋和質量評估

2）查看方式

# zcat查看gzip壓縮的文件 # head -n 8 顯示前8行文件內容（前8行代表2條序列）zcat filename.fq.gz | head -n 8# @SRR1039521.13952745/1 # TTCCTTCCTCCTCTCCCTCCCTCCCTCCTTTCTTTCTTCCTGTGGTTTTTTCCTCTCTTCTTC # + # HIJIIJHGHHIJIIIJJJJJJJJJJJJJJJJJJJJJIIJJFIDHIBGHJIHHHHHHFFFFFFE

3）格式說明：fastq文件每4行代表一條序列
第一行：記錄序列測序時所用儀器以及在測序通道中坐標信息，以@開頭；
第二行：測序的序列信息，以ATCGN表示，由于熒光信號干擾無法判斷是什么堿基時就用N表示；
第三行：通常一個+;
第四行：與第二行堿基信息一一對應，存儲測序堿基的質量值。

4）其他常用命令

# 計算read數 # wc -l: 計算行數 # bc -l: 計算器 (-l：浮點運算) # 為什么除以4，又除以1000000，計算的是million值echo "`zcat trt_N061011_1.fq.gz | wc-l` / (4*1000000)" | bc -l# 測序堿基數計算 zcat trt_N061011_1.fq.gz | awk'{if(FNR%4==0) base+=length}END{print base/10^9,"G";}'

awk的介紹見：常用和不太常用的awk命令

2.基因組FASTA文件

此文件可以從ensemble數據庫下載的（https://www.ensembl.org/info/data/ftp/index.html）， 一般選擇下載primary assemblyfasta（想知道為什么，點這里）。fasta文件用于序列存儲，可以是DNA或蛋白序列，在此FASTA文件存儲了基因組序列的信息。

序列名字行：以>符號開頭，記錄了該序列類型和所在基因組位置信息；

序列行（一行或多行）：序列信息，soft-masked基因組會把所有重復區和低復雜區的序列用小寫字母標出的基因組，小寫字母n表示未知堿基。

>1 dna_sm:chromosomechromosome:GRCh38:1:1:248956422:1 REF nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn ..... ttgggctggggcctggccatgtgtatttttttaaatttccactgatgattttgctgcatg gccggtgttgagaatgactgCGCAAATTTGCCGGATTTCCTTTGCTGTTCCTGCATGTAG TTTAAACGAGATTGCCAGCACCGGGTATCATTCACCATTTTTCTTTTCGTTAACTTGCCG ..... nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn# 通常要求序列名字行簡單為好，而且一般加chr作為開頭 # 給第一行添加chr標簽，并去掉其他多余信息 # 下面的寫法復雜了些，是為了避免給已經有chr信息的名字再加一次 # 幫助無腦操作 sed 's/^>\([^chr]\)/>chr\1/' ?Homo_sapiens.GRCh38.dna.primary_assembly.fa |cut -f 1 -d ' ' > GRCh38.fa

3. 基因組注釋文件gff和gtf

gff全稱General featureformat，主要是用來注釋基因組。gtf全稱Gene transfer format，主要是用來對基因進行注釋。兩者均是一個9列的基因信息注釋文件，前8列的信息幾乎一樣，區別在于第9列。具體可見歷史推文NGS基礎 - GTF/GFF文件格式解讀和轉換 在此不再贅述。

從ensemble下載的gtf文件前5行一般是以#開頭的注釋信息，后續分析中用不上需要去除，同時需要給第一列添加chr標簽（與基因組序列一致），可通過下面的命令對文件進行加工：

# grep 匹配查詢 -v 輸出不匹配的行 gunzip Homo_sapiens.GRCh38.94.gtf.gz -c |grep -v '^#' | sed '/^[^chr]/ s/^/chr/' >GRCh38.gtf

4. bed文件

分析過程中的bed文件一般代表區域信息，如表示Peak位置的bed文件，表示基因注釋的bed12文件。

• 表示基因注釋時，gtf/gff和bed文件的區別

1）gtf/gff文件一行表示一個exon/CDS等子區域，多行聯合表示一個gene；bed文件一行表示一個gene；
2）gtf文件中堿基位置定位方式是1-based，而bed中堿基定位方式是0-based，如下圖所示。

• bed文件每一行對應信息

• 
  

必須包含的3列信息：
1）chrom：染色體名字 (e.g.chr3, chrY, chr2_random或者scaffold10671)。
2）chromStart：基因在染色體或scaffold上的起始位置（0-based）。
3）chromEnd：基因在染色體或scaffold上的終止位置 （前閉后開）。

可選的9列信息：
4）name：bed文件的行名。
5）score：本條基因在注釋數據集文件中的評分（0-1000），在Genome Browser中會根據不同區段的評分顯示對應的陰影強度（評分越高灰度越高）。
6）strand：鏈的方向+、-或. (.表示不確定鏈的方向)
7）thickStart：CDS區（編碼區）的起始位置，即起始密碼子的位置。
8）thickEnd：The endingposition at which the feature is drawn thickly (for example the stop codon ingene displays).
9）itemRgb：RGB顏色值（如：255,0,0），方便在GenomeBrowser中查看。
10）blockCount：bed行中外顯子的數目。
11）blockSizes：逗號分割的列，數目與blockCount值對應，每個數表示對應外顯子的堿基數。
12）blockStarts：逗號分割的列，數目與blockCount值對應，每個數表示對應外顯子的起始位置（數值是相對ChromStart計算的）。

5. sam和bam文件

sam文件全稱The SequencingAlignment/Map Format，是Alignment/Map步驟bwa/STAR/HISAT2等軟件對結果的標準輸出文件，用于存儲reads比對到參考基因組的比對結果，是一個純文本格式，文件一般較大。為了節省硬盤存儲，一般使用其高效壓縮的二進制格式bam文件。

利用samtools view的-b參數就能把sam文件轉為bam文件。

1）sam文件查看方式
在linux終端直接用less即可進行查看；

2）bam文件查看方式
需要借助samtools view工具進行查看

samtools view filename.bam | less -S samtools view -h filename.bam | less -S

NGS分析中大多數文件都是由header和record兩部分組成，加上-h參數后可以將header顯示出來，默認是不顯示的。

@HD ? ?VN:1.5 ?SO:coordinate @SQ ? ?SN:chr1 LN:248956422 @SQ ? ?SN:chr10 ? ? ? ?LN:133797422 ...... @SQ ? ?SN:chrKI270392.1 ? ? ? ?LN:971 @SQ ? ?SN:chrKI270394.1 ? ? ? ?LN:970 @RG ? ?ID:BH_H3K27ac_2 LB:BH_H3K27ac_2 SM:BH_H3K27ac_2 @PG ? ?ID:bwa ?PN:bwa ?VN:0.7.15-r1140 CL:bwa mem -M -t 8 -R@RG\tID:BH_H3K27ac_2\tLB:BH_H3K27ac_2\tSM:BH_H3K27ac_2\tPL: /MP @PG ? ?ID:MarkDuplicates ? ? ?VN:1.138(aa51703435dc6a423013e74e56b0b68405facd79_1439324166) ? CL:picard.sam.markduplicates. K00141:244:HVL3NBBXX:8:2119:27235:3145399 ? ? ?chr1 ? ?10016 ?32 ? ? ?115M ? ?= ? ? ?10016 ? 115 ? ? CCCTAACCCTAACCCTAACCC K00141:244:HVL3NBBXX:8:2119:27235:31453147 ? ? chr1 ? ?10016 ?32 ? ? ?115M ? ?= ? ? ?10016 ? -115 ? CCCTTACCCTAACCCTAACCC

• header內容
  @HD：是必須的標準文件頭，包含版本信息；
  @SQ：參考序列染色體名字和長度信息 (SN:scaffold name; LN: length)；
  @RG：重要read group信息，通常包含測序平臺，測序文庫和樣本ID等信息，分析時用于區分不同樣本（重測序時用到）；
  @PG：生成此文件的操作過程和參數信息 （program）。

• record內容
  每一行就是一條read比對上參考基因組的信息，總共12列，用tab鍵分割。

# 1. read名稱； # 2. 比對信息位flag值； # 3. 參考序列染色體編號； # 4. 5′端起始位置； # 5. MAPQ：mapping quality，描述比對的質量，數字越大，特異性越高； # 6. CIGAR字符串，記錄插入、刪除、錯配等信息； # 7. 配對read所比對到的染色體，僅雙端測序的數據才有； # 8. 配對read所比對到的位置，僅雙端測序的數據才有； # 9. 插入片段的長度，僅雙端測序的數據才有； # 10. read序列； # 11. read質量值； # 12. 12列以后的信息都是metadata，程序用標記

sam文件中第二列flag信息很重要，下面做進一步解釋。

利用samtools flagstat工具可以查看bam文件中比對的flag信息，并輸出比對的統計結果。

samtools flagstat *.bam

flag一共有12個標簽，使用16進制數表示，每個標簽值是2^(n-1)，其中n<=12，每個值有其對應的唯一解釋含義，具體見下圖。

你會發現隨機挑選幾個值做加和運算，他們的結果都是唯一的，所以在bam文件中第二列flag的值代表這條序列符合下圖所示條件的值的和。

所以根據這個值我們可以判斷這條序列是雙端測序還是單端測序；如果是雙端測序，reads來自左端還是右端。比如65 只能是1和64組成，代表這個序列是雙端測序，而且是read1。

每次轉換很頭疼？別擔心，網上有很多解碼flag含義的在線工具，如SAM Format（網址：https://www.samformat.info/sam-format-flag）

輸入flag的值，解析工具會返回匹配上的信息。如果是單端測序，flag值都是偶數。

如果是雙端測序，工具會幫我們把另外一端序列的flag值返回，并且將這些數字情況大致分為5類，在右側進一步顯示這個值對應的含義。

6. wig、bigwig和bedgraph文件

上述bam和sam文件可以幫助我們探索reads在參考基因組中的比對情況，導入基因組瀏覽器查看比對狀態和突變信息。而wiggle(簡稱wig)、bigwig(簡寫bw)以及bedgraph(簡寫bdg)只包含區域和區域的覆蓋度信息，文件更小，用于可視化更方便，可以導入基因組瀏覽器（Genome Browser）中進行可視化，以查看reads在參考基因組各個區域的覆蓋度并檢測測序深度。這幾個文件在ChIP-seq數據分析Call Peak階段會生成，可以利用IGV等工具進行可視化，方便查看組蛋白修飾的程度。

• wiggle：展示區域的密度或強度信息，如GCpercent, probability scores, and transcriptome data.

variableStep chrom=chr2 300701 12.5 300702 12.5 300703 12.5 300704 12.5 300705 12.5fixedStep chrom=chr3 start=400601 step=100span=5 11 22 33

• bedGraph: bed與wig的結合，更省空間和靈活，展示信息與wig類似。(bedGraph的格式一般有四列，Chr、start、end和value，并且坐標是以0為起始左閉右開)

chromA chromStartA ?chromEndA ?dataValueA chromB chromStartB ?chromEndB ?dataValueB

• bigWig: wig文件的二進制壓縮格式，可通過wigToBigWig工具轉換

推薦大家閱讀UCSC官網對這幾個文件的詳細解釋：

• wiggle(WIG)：https://genome.ucsc.edu/goldenPath/help/wiggle.html

• bedGraph：https://genome.ucsc.edu/goldenPath/help/bedgraph.html

• bigWig：http://genome.ucsc.edu/goldenPath/help/bigWig.html

測試數據可視化參見歷史推文：

• 測試數據可視化（一）

• 測序數據可視化 （二）-IGV

• 測序數據可視化 （三） - UCSCgenomebrowser

• 測序數據可視化 （四）-Epigenomebrowser

• IGV基因組瀏覽器可視化高通量測序數據

• 本地安裝UCSC基因組瀏覽器

• 高通量數據分析必備|基因組瀏覽器使用介紹 -1

• 高通量數據分析必備|基因組瀏覽器使用介紹 -2

• 高通量數據分析必備|基因組瀏覽器使用介紹 -3

往期精品(點擊圖片直達文字對應教程)

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的生信分析过程中这些常见文件的格式以及查看方式你都知道吗？的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。