WES測(cè)序分析：fastq(bwa)bam(samtools+bcftools)vcf

RNA-Seq測(cè)序分析：fastq(hisat2)bam(FeatureCounts)counts(DEseq2)DEG

ChIP-Seq數(shù)據(jù)分析：fastq(bowtie2)bam(deeptools+macs2)bed(homer/meme)motif

1 fasta fastq

fasta序列：

1、以>開(kāi)頭

2、gi | gi號(hào) | 來(lái)源標(biāo)識(shí) | 序列標(biāo)識(shí)

3、序列描述

4、序列中允許空格，換行，空行，直到下一個(gè)>，表示該序列結(jié)束

image.png

查看序列數(shù)量
$ cat Arabidopsis\ arenosaGCA_9028502551.fna | grep '>' | wc
 121101 1332111 12336289

fastq序列：序列+質(zhì)量

1、@開(kāi)頭，唯一的序列ID標(biāo)識(shí)符，可選序列描述內(nèi)容，以空格分開(kāi)。

2、序列字符

3、+開(kāi)頭，空或加第一行@后相同內(nèi)容

4、堿基質(zhì)量字符，這一行字符數(shù)與第二行字符數(shù)相同。

image.png

先安裝fastqc軟件

進(jìn)入rnaseq環(huán)境
conda activate rnaseq

搜索fastqc軟件
(rnaseq) root 21:12:16 /home/kaoku/data/fastq
$ mamba search fastqc

安裝fastqc
$ mamba install fastqc

取前100行為臨時(shí)文件tmp.fq分析
$ cat hcc1395_normal_rep1_r1.fastq | head -100 > tmp.fq

fastqc軟件分析tmp.fq
$ fastqc tmp.fq
Started analysis of tmp.fq
Analysis complete for tmp.fq

會(huì)得到.html和.zip的分析結(jié)果
$ ls
hcc1395_normal_rep1_r1.fastq  tmp_fastqc.html  tmp_fastqc.zip  tmp.fq

進(jìn)一步查看tmp.fq
$ cat tmp.fq | paste - - - - | less -S

image.png

分析的結(jié)果：

image.png

$ cat tmp.fq | paste - - - - | wc
  25000  100000 8688933

查看第二列的首個(gè)堿基分布情況
$ cat tmp.fq | paste - - - - | cut -f 2 | cut -c1 | sort | uniq -c
   2223 A
   7905 C
  13095 G
      8 N
   1769 T

第二位堿基情況
$ cat tmp.fq | paste - - - - | cut -f 2 | cut -c2 | sort | uniq -c
   2931 A
   6630 C
   7026 G
   8413 T

為了讓結(jié)果更可靠增加了堿基的數(shù)目
$ cat hcc1395_normal_rep1_r1.fastq | head -100000 > tmp.fq

2 SAM BAM

SAM：

1、標(biāo)頭注釋文件 以@開(kāi)頭：@HD：符合標(biāo)準(zhǔn)的版本、對(duì)比序列，@SQ：參考序列說(shuō)明，@PG：使用的對(duì)比程序說(shuō)明，LN：參考序列的長(zhǎng)度

2、對(duì)比結(jié)果部分

現(xiàn)成的SAM文件：

$ less -S tumor_rep1.sam

image.png

我們把先前的tmp.fq文件比對(duì)成SAM文件。

從安裝軟件開(kāi)始

安裝bwa
$ mamba install bwa
或者
$ conda install -c bioconda bwa

查看bwa目錄
$ bwa
Command: index         index sequences in the FASTA format
         mem           BWA-MEM algorithm
         fastmap       identify super-maximal exact matches
         pemerge       merge overlapping paired ends (EXPERIMENTAL)
         aln           gapped/ungapped alignment
         samse         generate alignment (single ended)
         sampe         generate alignment (paired ended)
         bwasw         BWA-SW for long queries

         shm           manage indices in shared memory
         fa2pac        convert FASTA to PAC format
         pac2bwt       generate BWT from PAC
         pac2bwtgen    alternative algorithm for generating BWT
         bwtupdate     update .bwt to the new format
         bwt2sa        generate SA from BWT and Occ

采用mem子菜單，開(kāi)始比對(duì)（下載參考基因組后需要先構(gòu)建索引）
$ bwa men hg38 tmp.fq > tmp.sam

比對(duì)結(jié)束后顯示

SAM文件轉(zhuǎn)BAM文件（就是壓縮）
$ samtools view -bS tmp.sam > tmp.bam 即可

查看BAM文件
$ samtools view tmp.bam

3 gff gtf

gtf 空格隔開(kāi) 9列，參考序列ID、注釋來(lái)源、類(lèi)型、開(kāi)始位點(diǎn)、結(jié)束位點(diǎn)、得分、正負(fù)鏈、步進(jìn)、屬性。

gff =號(hào)隔開(kāi)

4 Bigwig Wiggle

測(cè)序深度、覆蓋度

5 bed

略

6 vcf

查看變異，測(cè)序深度等

samtools mpileup -ugf hg38.fa tmp.bam

samtools mpileup -ugf hg38.fa tmp.bam | bcftools call -vm0 z -o

7 Linux習(xí)題

7.1 fa-fq習(xí)題

$ ls
combined_reads.bam  longreads.fq  reads_1.fq  reads_2.fq  simulate.pl

1）統(tǒng)計(jì)reads_1.fq 文件中共有多少條序列信息
$ wc -l reads_1.fq
40000 reads_1.fq

$ cat reads_1.fq | paste - - - - | wc
  10000   40000 2285692

2）輸出所有的reads_1.fq文件中的標(biāo)識(shí)符(即以@開(kāi)頭的那一行)
$ cat reads_1.fq | paste - - - - | cut -f1
或 $ awk '{if(NR%1==0)print}' reads_1.fq （NR行號(hào)）

3）輸出reads_1.fq文件中的 所有序列信息(即每個(gè)序列的第二行)
$ cat reads_1.fq | paste - - - - | cut -f2

4）輸出以‘+’及其后面的描述信息(即每個(gè)序列的第三行)
$ cat reads_1.fq | paste - - - - | cut -f3

5）輸出質(zhì)量值信息(即每個(gè)序列的第四行)
$ cat reads_1.fq | paste - - - - | cut -f4
或 $ awk '{if(NR%4==0)print}' reads_1.fq

6）計(jì)算reads_1.fq 文件含有N堿基的reads個(gè)數(shù)
$ awk '{if(NR%4==2)print}' reads_1.fq | grep N | wc
   6429    6429  782897
或$ awk '{if(NR%4==2)print}' reads_1.fq | grep -c N

7）統(tǒng)計(jì)文件中reads_1.fq文件里面的序列的堿基總數(shù)
$ awk '{if(NR%4==2)print}' reads_1.fq | grep -o [ATCGN] | wc
1088399 1088399 2176798
或者把所有行數(shù)相加
$ awk '{if(NR%4==2)print length}' reads_1.fq | paste -s -d + | bc
1088399

8）計(jì)算reads_1.fq 所有的reads中N堿基的總數(shù)
$ awk '{if(NR%4==2)print}' reads_1.fq | grep -o [N] | wc
  26001   26001   52002

9）統(tǒng)計(jì)reads_1.fq 中測(cè)序堿基質(zhì)量值恰好為Q20的個(gè)數(shù)
$ awk '{if(NR%4==0)print}' reads_1.fq | grep -o '5' | wc
  21369   21369   42738

10）統(tǒng)計(jì)reads_1.fq 中測(cè)序堿基質(zhì)量值恰好為Q30的個(gè)數(shù)
$ awk '{if(NR%4==0)print}' reads_1.fq | grep -o '?' | wc
  21574   21574   43148

11）統(tǒng)計(jì)reads_1.fq 中所有序列的第一位堿基的ATCGN分布情況
$ awk '{if(NR%4==2)print}' reads_1.fq | cut -c1 | sort | uniq -c
   2184 A
   2203 C
   2219 G
   1141 N
   2253 T

12）將reads_1.fq 轉(zhuǎn)為reads_1.fa文件(即將fastq轉(zhuǎn)化為fasta)
$ cat reads_1.fq | paste - - - - | cut -f1,2 | tr '\t' '\n'| tr '@' '>' > reads_1.fa

13）統(tǒng)計(jì)上述reads_1.fa文件中共有多少條序列
$ wc reads_1.fa
  20000   20000 1167293 reads_1.fa

14）計(jì)算reads_1.fa文件中總的堿基序列的GC數(shù)量
$ grep -o [GC] reads_1.fa | wc
 529983  529983 1059966

15）刪除 reads_1.fa文件中的每條序列的N堿基
$ cat reads_1.fa | tr -d 'N'

16）刪除 reads_1.fa文件中的含有N堿基的序列
$ cat reads_1.fa | paste - - | grep -v N | tr '\t' '\n'

17）刪除 reads_1.fa文件中的短于65bp的序列
$ cat reads_1.fa | paste - - | awk '{if(length($2)>65)print}'| wc
   7076   14152  992399

18）刪除 reads_1.fa文件每條序列的前后五個(gè)堿基

19）刪除 reads_1.fa文件中的長(zhǎng)于125bp的序列

20）查看reads_1.fq 中每條序列的第一位堿基的質(zhì)量值的平均值

7.2 sam-bam習(xí)題

image.png

獲取文件
$  bowtie2 -x /home/kaoku/data/tmp/example/index/lambda_virus -1 reads_1.fq -2 reads_2.fq > tmp.sam                 

$ samtools view -bS tmp.sam > tmp.bam

1）統(tǒng)計(jì)共多少條reads(pair-end reads這里算一條)參與了比對(duì)參考基因組
$ cat tmp.sam | grep -v "^@" | wc
  20000  391929 7048498
答案是10000條

2）統(tǒng)計(jì)共有多少種比對(duì)的類(lèi)型(即第二列數(shù)值有多少種)及其分布。
$ cat tmp.sam | grep -v "^@" | cut -f2 | sort | uniq -c | sort -k1,1nr
   4650 163
   4650 83
   4516 147
   4516 99
    213 141
    213 77
    165 137
    165 69
    153 133
    153 73
    136 165
    136 89
    125 101
    125 153
     24 129
     24 65
     16 113
     16 177
      2 161
      2 81
搜一下163什么意思。

3）篩選出比對(duì)失敗的reads，看看序列特征。
$ cat tmp.sam | grep -v "^@" | awk '{if($6=="*")print}'|wc
   1005   12608  255140

查看序列$ cat tmp.sam | grep -v "^@" | awk '{if($6=="*")print$10}'
有很多N

4）比對(duì)失敗的reads區(qū)分成單端失敗和雙端失敗情況，并且拿到序列ID
雙端沒(méi)比對(duì)上$ cat tmp.sam | grep -v "^@" | awk '{if($6=="*")print$1}'| sort | uniq -c | grep -w 2
單端沒(méi)比對(duì)上$ cat tmp.sam | grep -v "^@" | awk '{if($6=="*")print$1}'| sort | uniq -c | grep -w 1

5）篩選出比對(duì)質(zhì)量值大于30的情況（看第5列）
$ cat tmp.sam | grep -v "^@" | awk '{if($5>30)print$0}' | less -S

6）篩選出比對(duì)成功，但是并不是完全匹配的序列
$ cat tmp.sam | grep -v "^@" | awk '{if($6!="*")print$6}' | grep "[IDNXSHP]"
D居多

7）篩選出inset size長(zhǎng)度大于1250bp的 pair-end reads
$ cat tmp.sam | grep -v "^@" | awk '{if($7>1250)print}'| less -S

8）統(tǒng)計(jì)參考基因組上面各條染色體的成功比對(duì)reads數(shù)量
$ cut -f 3 tmp.sam | sort -u
*
gi|9626243|ref|NC_001416.1|
LN:48502
PN:bowtie2
SO:unsorted
只有一條

9）篩選出原始fq序列里面有N的比對(duì)情況
$ awk '{if($10~"N")print}' tmp.sam | less -s

10）篩選出原始fq序列里面有N，但是比對(duì)的時(shí)候卻是完全匹配的情況
$ awk '{if($10~"N")print}' tmp.sam | awk '{if($6!~"[NIDSHP]")print}'| grep -v * | less -S

11）sam文件里面的頭文件行數(shù)
$ head sam.tmp
$ grep @ |wc

12）sam文件里每一行的tags個(gè)數(shù)一樣嗎
$ cut -f 12-1000 tmp.sam
不一樣
13）sam文件里每一行的tags個(gè)數(shù)分別是多少個(gè)

14）sam文件里記錄的參考基因組染色體長(zhǎng)度分別是？
頭文件有記錄
$ head tmp.sam
@HD     VN:1.5  SO:unsorted     GO:query
@SQ     SN:gi|9626243|ref|NC_001416.1|  LN:48502

15）找到比對(duì)情況有insertion情況的
$ awk '{if($6~"I")print}' tmp.sam | less -S

16）找到比對(duì)情況有deletion情況的
$ awk '{if($6~"D")print}' tmp.sam | less -S

17）取出位于參考基因組某區(qū)域的比對(duì)記錄，比如 5013到50130 區(qū)域
$ awk '{if($4>5013 && $4<50130)print}' tmp.sam | less -S

18）把sam文件按照染色體以及起始坐標(biāo)排序
sort

19）找到 102M3D11M 的比對(duì)情況，計(jì)算其reads片段長(zhǎng)度。
$ grep 102M3D11M tmp.sam
r284    83      gi|9626243|ref|NC_001416.1|     20953   42      102M3D11M       =       20869   -200    CGCGTCACCTGACGCACTGAATACGCTGAATGAACTGGCCGCAGCGCTCGGGAATGATCCAGATTTTGCTACCACCATGACTAACGCGNTTGCGGGTAAACAGAAGAATGCGA    6C1-F-'BB,0%-")!#?4,'?;A=B76$"B%3B!3/#%!36H*;D>/HF78>2720''EDD!H1E%?/8D2G/(D4+@"+H8D2EA6EEA:';:/<5;/HH:<.5&!)"E!.    AS:i:-15        XN:i:0  XM:i:1  XO:i:1  XG:i:3  NM:i:4  MD:Z:88C13^ACC11        YS:i:-3 YT:Z:CP

reads:CGCGTCACCTGACGCACTGAATACGCTGAATGAACTGGCCGCAGCGCTCGGGAATGATCCAGATTTTGCTACCACCATGACTAACGCGNTTGCGGGTAAACAGAAGAATGCGA

20）安裝samtools軟件后使用samtools軟件的各個(gè)功能?chē)L試把上述題目重新做一遍。

下一篇我們將補(bǔ)充Linux基礎(chǔ)習(xí)題20問(wèn)。

我們下一篇再見(jiàn)！