刘小泽写于2018年8月，发送于19.7.5
基因组的那些事儿--基础：https://www.jianshu.com/p/bf871522ea20
基因组的那些事儿(二)：https://www.jianshu.com/p/58a895a81ec6
基因组的那些事儿(三)-准备工作：https://www.jianshu.com/p/29ed65d3608b
4 质量控制(放在kpgp_wgs/qc下)
使用fastp进行质控，速度十分快，我之前测试了24G的双端测序文件，不到六分钟处理完成。下载的KPGP原始数据共68G，6个双端测序文件，平均一个用时1300秒
因为是从数据库中下载的现成的测序文件，质量已经调到了最好，用fastp质控后，改动不是很大【如果是测序质量差的文件，在fastp过滤后会有明显差别】如图是我们的KPGP的原始数据
可以自己统计GC、Q20、Q30
#!/usr/bin/perl
opendir (DIR, "./") or die "can't open the directory!";
@dir = readdir DIR;
foreach $file ( sort @dir)
{
next unless -d $file;
next if $file eq '.';
next if $file eq '..';
$total_reads= `grep '^Total' ./$file/fastqc_data.txt`;
$total_reads=(split(/\s+/,$total_reads))[2];
$GC= `grep '%GC' ./$file/fastqc_data.txt`;
$GC=(split(/\s+/,$GC))[1];
chomp $GC;
open FH , "<./$file/fastqc_data.txt";
while (<FH>)
{
next unless /#Quality/;
while (<FH>)
{
@F=split;
$hash{$F[0]}=$F[1];
last if />>END_MODULE/;
}
}
$all=0;$Q20=0;$Q30=0;
$all+=$hash{$_} foreach keys %hash;
$Q20+=$hash{$_} foreach 0..20;
$Q30+=$hash{$_} foreach 0..30;
$Q20=1-$Q20/$all;
$Q30=1-$Q30/$all;
print "$file\t$total_reads\t$GC\t$Q20\t$Q30\n";
#print "$all\n";
#
}
5 比对：终其一生，只为遇见你(放在kpgp_wgs/align下)
何为终其一生？
每条测序reads长度是150bp，只有reads这些ATCG字母，对我们是没有用的，起码要知道这些ATCG的排列在基因组的什么地方吧。人类的参考基因组30亿个碱基，如果要150一段150一段去观察的话，真的，一条reads要终其一生，只为遇见那个和他匹配的那个她。更何况，一个全基因组测序下来就有至少6亿条reads😱
但！我们有不怕苦不怕累的计算机！可以使用软件bwa或者bowtie，将fastq测序文件与参考基因组做对比，得到的结果是sam格式的文件【如下图】，其中可以看到每条reads在参考基因组的位置，在哪一条染色体，或者在这条染色体的哪个位置
5.1 Fastq to sam
为了得到这个结果，我们可以使用bwa软件进行操作
#格式：bwa mem -t线程数 -M处理同一个reads比对到参考基因组上不同位置的情况 index名称（要定位到前缀，不能只定位到文件夹名，比如这里倒数第二个hg19是文件夹名，最后一个hg19才是前缀）输入fastq1 输入fastq2 1>标准输出（sam格式）2>标准错误输出
for i in $(seq 1 6);do
bwa mem -t 10 -M ~/reference/index/bwa/hg19/hg19 ../raw_data/KPGP-00001_L${i}_R1.fq.gz ../raw_data/KPGP-00001_L${i}_R2.fq.gz 1>KPGP-00001_L${i}.sam 2>KPGP-00001_L${i}.bwa.align.log
done
【使用60线程运行，平均50分钟比对完一个】得到的结果是sam格式的，sam的全称是sequence alignment/map format，往往文件很大，为了节省空间又保证同样效果，他的孪生兄弟bam诞生【bam是sam的二进制文件，占用空间要小很多】
5.2 sam to bam
需要用samtools进行格式转换【转换依旧使用60线程，大概半个小时就处理好了，可以看出bam比sam节省了三分之二的空间】
#sam2bam -b设置输出为bam格式 -S 指定输入为sam
samtools view -b -S file.sam > file.bam
#bam2sam -h意思是加上header
samtools view -h file.bam > file.sam
5.3 sam/bam文件格式
有了文件，以后的大部分分析都是基于sam格式进行的，因此理解格式至关重要
记得这张图当时做了半个多小时
第3和第7列，可以用来判断某条reads是否成功比对到了基因组的染色体上，另外两条PE reads是否比对到同一条染色体；
第1，10，11列可以还原成测序数据的fastq格式
5.4 对bam文件进行排序=>sorted_bam
for i in `seq 1 6`;do
samtools sort -@ 15 -o KPGP-00001_L${i}.sorted.bam KPGP-00001_L${i}.bam
done
每个bam文件按照最左侧比对位置的先后进行排序，方便接下来的构建索引
5.5 对sorted_bam文件构建索引=>sorted_bam.bai
for i in `seq 1 6`;do
samtools index KPGP-00001_L${i}.sorted.bam
done
bam文件进行默认情况下的坐标排序后,才能进行index。目的是为了更快地访问bam文件，结果会生成.bai格式的索引文件
samtools index是针对bam用的，如果要对sam文件构建索引，要用samtools tabix命令
5.6 查看参考基因组每个位点的比对情况=>mpileup.txt
利用samtools的mpileup功能可以生成bcf文件，之后可以结合bcftools进行SNP和InDel的分析；
参数：-f：输入有索引的参考基因组fasta；
-g：输出二进制的bcf格式（不加-g就只生成文本文件，统计了针对不同的测序比对文件，参考序列中每个碱基位点的情况，每一行表示参考序列中某个碱基的比对结果）
for i in `seq 1 6`;do
samtools mpileup -f ~/reference/genome/hg19/hg19.fa KPGP-00001_L${i}.sorted.bam > KPGP-00001_L${i}.mpileup.txt
done
查看各个测序样本bam比对结果
for i in `seq 1 6`;do
samtools flagstat KPGP-00001_L${i}.sorted.bam >KPGP-00001_L${i}.flagstat.txt
done
以上几步操作每个测序比对文件会形成
5.7 将各个bam文件融合=>merge.bam
原来只是各个测序文件比对的结果，他们虽然都sort和index了，也有了各自与参考基因组的位点比对情况。但是我们不能从整体角度去查看，要想统一原本的“各自为政”，我们需要进行merge
samtools merge -@ 20 KPGP-00001.merge.bam *.sorted.bam
5.8 对merge后的bam也进行排序=>sorted.merge.bam
##sort
samtools sort -@ 20 -O bam -o KPGP-00001.sorted.merge.bam KPGP-00001.merge.bam
##index
samtools index KPGP-00001.sorted.merge.bam
5.9 对merge_bam进行拆分=>n个小bam
结果得到的
KPGP-00001.sorted.merge.bam是61G，这样的文件是不可能直接导入IGV查看的，有没有什么方法可以缩小文件，并按照染色体进行查看呢？
查看一个大文件不靠谱，那么可不可以将大文件拆分开来，拆成一个个的小文件呢？利用bamtools split是可以的。另外，如果程序随意给我们排序好的bam文件给拆了，那么之前的排序工作也没有任何意义，因此，我们可以指定参数-reference，让程序根据参考序列（也就是根据bam第三列的染色体编号） 进行拆分
bamtools split -in KPGP-00001.sorted.merge.bam -reference
下面这张图就是合并bam、合并后排序的bam以及拆分后的bam全家福
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
网友评论