WES测序分析：fastq(bwa)bam(samtools+bcftools)vcf

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

ChIP-Seq数据分析：fastq(bowtie2)bam(deeptools+macs2)bed(homer/meme)motif

1 fasta fastq

fasta序列：

1、以>开头

2、gi | gi号 | 来源标识 | 序列标识

3、序列描述

4、序列中允许空格，换行，空行，直到下一个>，表示该序列结束

image.png

查看序列数量
$ cat Arabidopsis\ arenosaGCA_9028502551.fna | grep '>' | wc
 121101 1332111 12336289

fastq序列：序列+质量

1、@开头，唯一的序列ID标识符，可选序列描述内容，以空格分开。

2、序列字符

3、+开头，空或加第一行@后相同内容

4、碱基质量字符，这一行字符数与第二行字符数相同。

image.png

先安装fastqc软件

进入rnaseq环境
conda activate rnaseq

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

安装fastqc
$ mamba install fastqc

取前100行为临时文件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

会得到.html和.zip的分析结果
$ ls
hcc1395_normal_rep1_r1.fastq  tmp_fastqc.html  tmp_fastqc.zip  tmp.fq

进一步查看tmp.fq
$ cat tmp.fq | paste - - - - | less -S

image.png

分析的结果：

image.png

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

查看第二列的首个碱基分布情况
$ 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

为了让结果更可靠增加了碱基的数目
$ cat hcc1395_normal_rep1_r1.fastq | head -100000 > tmp.fq

2 SAM BAM

SAM：

1、标头注释文件 以@开头：@HD：符合标准的版本、对比序列，@SQ：参考序列说明，@PG：使用的对比程序说明，LN：参考序列的长度

2、对比结果部分

现成的SAM文件：

$ less -S tumor_rep1.sam

image.png

我们把先前的tmp.fq文件比对成SAM文件。

从安装软件开始

安装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子菜单，开始比对（下载参考基因组后需要先构建索引）
$ bwa men hg38 tmp.fq > tmp.sam

比对结束后显示

SAM文件转BAM文件（就是压缩）
$ samtools view -bS tmp.sam > tmp.bam 即可

查看BAM文件
$ samtools view tmp.bam

3 gff gtf

gtf 空格隔开 9列，参考序列ID、注释来源、类型、开始位点、结束位点、得分、正负链、步进、属性。

gff =号隔开

4 Bigwig Wiggle

测序深度、覆盖度

5 bed

略

6 vcf

查看变异，测序深度等

samtools mpileup -ugf hg38.fa tmp.bam

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

7 Linux习题

7.1 fa-fq习题

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

1）统计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文件中的标识符(即以@开头的那一行)
$ cat reads_1.fq | paste - - - - | cut -f1
或 $ awk '{if(NR%1==0)print}' reads_1.fq （NR行号）

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

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

5）输出质量值信息(即每个序列的第四行)
$ cat reads_1.fq | paste - - - - | cut -f4
或 $ awk '{if(NR%4==0)print}' reads_1.fq

6）计算reads_1.fq 文件含有N碱基的reads个数
$ 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）统计文件中reads_1.fq文件里面的序列的碱基总数
$ awk '{if(NR%4==2)print}' reads_1.fq | grep -o [ATCGN] | wc
1088399 1088399 2176798
或者把所有行数相加
$ awk '{if(NR%4==2)print length}' reads_1.fq | paste -s -d + | bc
1088399

8）计算reads_1.fq 所有的reads中N碱基的总数
$ awk '{if(NR%4==2)print}' reads_1.fq | grep -o [N] | wc
  26001   26001   52002

9）统计reads_1.fq 中测序碱基质量值恰好为Q20的个数
$ awk '{if(NR%4==0)print}' reads_1.fq | grep -o '5' | wc
  21369   21369   42738

10）统计reads_1.fq 中测序碱基质量值恰好为Q30的个数
$ awk '{if(NR%4==0)print}' reads_1.fq | grep -o '?' | wc
  21574   21574   43148

11）统计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 转为reads_1.fa文件(即将fastq转化为fasta)
$ cat reads_1.fq | paste - - - - | cut -f1,2 | tr '\t' '\n'| tr '@' '>' > reads_1.fa

13）统计上述reads_1.fa文件中共有多少条序列
$ wc reads_1.fa
  20000   20000 1167293 reads_1.fa

14）计算reads_1.fa文件中总的碱基序列的GC数量
$ 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文件每条序列的前后五个碱基

19）删除 reads_1.fa文件中的长于125bp的序列

20）查看reads_1.fq 中每条序列的第一位碱基的质量值的平均值

7.2 sam-bam习题

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）统计共多少条reads(pair-end reads这里算一条)参与了比对参考基因组
$ cat tmp.sam | grep -v "^@" | wc
  20000  391929 7048498
答案是10000条

2）统计共有多少种比对的类型(即第二列数值有多少种)及其分布。
$ 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）筛选出比对失败的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）比对失败的reads区分成单端失败和双端失败情况，并且拿到序列ID
双端没比对上$ cat tmp.sam | grep -v "^@" | awk '{if($6=="*")print$1}'| sort | uniq -c | grep -w 2
单端没比对上$ cat tmp.sam | grep -v "^@" | awk '{if($6=="*")print$1}'| sort | uniq -c | grep -w 1

5）筛选出比对质量值大于30的情况（看第5列）
$ cat tmp.sam | grep -v "^@" | awk '{if($5>30)print$0}' | less -S

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

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

8）统计参考基因组上面各条染色体的成功比对reads数量
$ cut -f 3 tmp.sam | sort -u
*
gi|9626243|ref|NC_001416.1|
LN:48502
PN:bowtie2
SO:unsorted
只有一条

9）筛选出原始fq序列里面有N的比对情况
$ awk '{if($10~"N")print}' tmp.sam | less -s

10）筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况
$ awk '{if($10~"N")print}' tmp.sam | awk '{if($6!~"[NIDSHP]")print}'| grep -v * | less -S

11）sam文件里面的头文件行数
$ head sam.tmp
$ grep @ |wc

12）sam文件里每一行的tags个数一样吗
$ cut -f 12-1000 tmp.sam
不一样
13）sam文件里每一行的tags个数分别是多少个

14）sam文件里记录的参考基因组染色体长度分别是？
头文件有记录
$ head tmp.sam
@HD     VN:1.5  SO:unsorted     GO:query
@SQ     SN:gi|9626243|ref|NC_001416.1|  LN:48502

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

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

17）取出位于参考基因组某区域的比对记录，比如 5013到50130 区域
$ awk '{if($4>5013 && $4<50130)print}' tmp.sam | less -S

18）把sam文件按照染色体以及起始坐标排序
sort

19）找到 102M3D11M 的比对情况，计算其reads片段长度。
$ 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软件的各个功能尝试把上述题目重新做一遍。

下一篇我们将补充Linux基础习题20问。

我们下一篇再见！