Skip to content

Latest commit

 

History

History
221 lines (164 loc) · 7.34 KB

测序数据评估.md

File metadata and controls

221 lines (164 loc) · 7.34 KB
Raw

目录

测序数据评估

测序数据评估

测序平台简介 目录

HiSeq X Ten 目录

HiSeq X Ten是Illumina于2014年推出的最新测序系统，其功能定位为工厂规模的测序系统，实现了Illumina测序仪迄今为止最高的测序通量和最低的测序成本。HiSeq X Ten系统由10台超高通量测序仪HiSeq X组成，测序读长为2×150bp，单台仪器每次运行可产出高达1.8Tb的数据，运行时间在三天以内，即每台仪器每天产出约600Gb的数据。10台仪器同时运行时，每周至少可完成320个人类基因组测序（以30×覆盖度计算），每年完成的数量可超过18000个。因此，HiSeq X Ten将使研究人员更易于开展大规模人类基因组测序，并将有利于深入挖掘与癌症及复杂疾病相关的遗传变异。

md5校验 目录

检查数据传输的完整性就是md5校验，看看数据在拷贝过程中有没有意外的损坏。

一般传输数据之前，会用md5命令来生成各个文件的md5值，并将md5值保存到MD5.txt文件中，然后传输数据之后，需要自行用md5sum -c MD5.txt 来校验文件里面记录的文件的完整性，如果显示都是OK，说明文件拷贝传输过程是没有问题的！但这个过程会耗费大量的磁盘读写，磁盘读写能力是有限的，所以开多个进程并不能加快这一过程。

注意：MD5.txt文件必须与对应的源文件保存在同一路径下

数据组织形式 目录

数据组织形式，即如何将原始数据与分析过程中产生的中间数据，按照性质分门别类地保存在对应的文件夹下，以方便数据的获取与管理

建议按照以下文件夹结构组织数据：

Cleandata
	|----Data	# fastq测序数据
	|----QC		# FastQC输出的结果文件
	|----Map	# Mapping过程的输出
	|----Stat	# 饱和度分析过程输出的饱和度分析统计结果

测序数据质量控制 目录

使用fastqc进行质控

$ fastqc  -o output dir [-(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
  • -o 是用来指定输出文件的目录，注意是这里是不能自动新建目录的
  • -noextract 输出的结果是.zip文件，默认自动解压缩，-noextract则不解压缩
  • -f 用来强制指定输入文件格式，默认会自动检测
  • -c 用来指定一个contaminant文件，fastqc会把overrepresented sequences往这个contaminant文件里搜索

所有的fastq.gz文件用fastqc软件处理得到的测序质量检测报告是一个html文件加上一个文件夹（若未自动解压则为一个html文件加上一个fastqc.zip文件）

批量解压

$ ls *zip|while read id;do unzip $id;done

统计fastqc结果 目录

用fastqc_stat.pl脚本统计fastqc结果。要提取的数据来自于解压后文件夹中的fastqc_data.txt文件

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
$total_reads=  `grep '^Total' ./$file/fastqc_data.txt`;
$total_reads=(split(/\s+/,$total_reads))[2];
# 提取%GC
$GC= `grep '%GC' ./$file/fastqc_data.txt`;
$GC=(split(/\s+/,$GC))[1];
chomp $GC;

# 提取Q20，Q30
## 读入Per sequence quality scores部分的信息，保存成哈希
open FH , "<./$file/fastqc_data.txt";
while (<FH>)
    {
    next unless /#Quality/;
    while (<FH>)
        {
        @F=split;
        $hash{$F[0]}=$F[1];
        last if />>END_MODULE/;
        }
    }
## 统计Q20，Q30
$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";
}

在FastQC的结果输出文件夹中运行该脚本

$ perl fastqc_stat.pl

饱和度分析 目录

基于测序深度 目录

需要先下好RefSeq作为参考序列，从UCSC上下载，下载地址：

http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/refMrna.fa.gz

用Bowtie2做好索引

$ nohup bowtie2-build --threads refMrna.fa refMrna.fa >refMrna.bwt_index.log 2>&1

将 reads 比对到 RefSeq 上

wd='..'
dir_data='..'
ls $dir_data/Path/To/Data | while read i
do
	{
	read1=`ls $dir_data/Path/To/Data/$i/*1.clean.fq.gz`
	read2=`ls $dir_data/Path/To/Data/$i/*2.clean.fq.gz`
	sampleName=`basename $read1 _1.clean.fq.gz`
	bowtie2 -p 8 --local -x $wd/Ref/hg19/RefSeq/refMrna.fa -1 $read1 -2 $read2 2>$dir_data/Path/To/Map/RefSeq_${sampleName}_map.stat | samtools sort -@ 8 -O bam -o $dir_data/Path/To/Map/RefSeq_${sampleName}.sort.bam 2>$dir_data/Path/To/Map/RefSeq_${sampleName}_map.log
	} &
done
wait

用Saturation_Depth_Stat.sh脚本进行饱和度分析（依据测序深度）

先准备好TotalReads.txt文本，共两列，一列为样本名，另一列为对应样本的reads数（双端总和）

NP006_RRS05401	104681638
NP007_RRS05402	123032856
OMI003_RRS05395	114242868
OMI024_RRS05396	119350468
STA002_RRS05397	108758530
STA003_RRS05398	108015800
STA027_RRS05399	104833546
STA031_RRS05400	112848138

Saturation_Depth_Stat.sh脚本：

dir_data='..'

# 判断统计结果输出的文件夹是否已经创建
if [ ! -d $dir_data/Path/To/Stat ]
then
	mkdir $dir_data/Path/To/Stat
fi

cat $dir_data/Path/To/TotalReads.txt | while read i
do
	{
	Sample=`echo $i | awk '{print $1}'`
	Total_reads=`echo $i | awk '{print $2}'`
	export Total_reads
	Exome_length=50000000
	Reads_length=150
	Depth1_reads=$[$Exome_length/$Reads_length]
	export Depth1_reads
	Max_depth=$[$Total_reads*$Reads_length/$Exome_length]
	if [ -f $dir_data/Path/To/Stat/${Sample}_depth.txt ]
	then
		rm $dir_data/Path/To/Stat/${Sample}_depth.txt
	fi
	for ((j=5;j<=Max_depth;j+=5))
	do
		Current_depth=$j
		export Current_depth
		echo -ne "$Current_depth\t" >>$dir_data/Path/To/Stat/${Sample}_depth.txt
		samtools view $dir_data/Path/To/Map/RefSeq_${Sample}.sort.bam | perl -ne 'chomp;next if (/^\@/);if (rand()<($ENV{"Current_depth"}*$ENV{"Depth1_reads"}/$ENV{"Total_reads"})){print "$_\n";}' | cut -f 3 | sort | uniq | wc -l >>$dir_data/Path/To/Stat/${Sample}_depth.txt
	done
} &
done
wait

批量执行饱和度分析

$nohup bash Saturation_Depth_Stat.sh &

参考资料：

(1) HiSeq X Ten

(2) 生信技能树《直播我的基因组》系列文章

(3) 写脚本对fastqc的结果进行统计咯！