|
其实我最害怕的就是写一个系列的文章，就像大学上课一样，刚开始的时候人来的很多，到最后人越来越少。 前情提要： 这次的目的是：
数据解压之前下载了所有的数据，但只有样本9~15才是mRNA-Seq测序结果，其中9-11为人类293个AKAP9敲除细胞，12-15为小鼠的AKAP9敲除细胞。也就是只要解压9~15就行了，即是SRR3589956 ~ SRR3589962.
所以基本上命令即使如下用法, 如果你觉得自己的空间够大就不需要用到 for id in `seq 56 62`do fastq-dump --gzip --split-3 -O /mnt/f/Data/RNA-Seq -A SRR35899${id}done压缩的文件不要着急解压，有很多bash命令能够直接用于压缩文件，如 zcat SRR3589956_1.fastq.gz | head -n 4用这些z-tools能够节省大量磁盘空间。 QC basic concept高通量测序之所以能够能够达到如此高的通量的原因就是他把原来几十M，几百M，甚至几个G的基因组通过物理或化学的方式打算成几百bp的短序列，然后同时测序。 因为测序过程是边合成边测序（SBS），所以在建库的时候，短序列两段会加一些固定的碱基用于桥式PCR扩增，这些固定的碱基就是adapter(接头)。一般而言，还可以在接头加一些tag（index），用于标识这个read来自于哪个物种。目前的单细胞测序为了省钱，譬如10X genomic技术，都是在一个pool里面加多种接头。 在测序过程中，机器会对每次读取的结果赋予一个值，用于表明它有多大把握结果是对的。从理论上都是前面质量好，后面质量差。并且在某些GC比例高的区域，测序质量会大幅度降低。 因此，我们在正式的数据分析之前需要对分析结果进行质控。Jimmy大神就发帖专门指出”要充分了解你的测序数据—论QC的重要性'，http://www./thread-324-1-1.html 。 Fastq文件格式说明FASTQ文件每个序列通常为4行，分别为：
FASTQ的文件示例： @DJB775P1:248:D0MDGACXX:7:1202:12362:49613 1:Y:18:ATCACGTGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA+JJJJJIIJJJJJJHIHHHGHFFFFFFCEEEEEDBD?DDDDDDBDDDABDDCA第一行序列名称
在v1.4之前
第三行质量序列格式
不同版本的碱基质量Q和碱基错误率P的关系如下 FastQC质量报告质量控制的软件很多，但是目前主要以fastqc为主。常见的用法： fastqc seqfile1 seqfile2 .. seqfileN常用参数：-o： 输出路径--extract: 输出文件是否需要自动解压 默认是--noextract-t: 线程， 和电脑配置有关，每个线程需要250MB的内存-c: 测序中可能会有污染， 比如说混入其他物种-a: 接头-q: 安静模式FastQC有两种方式分析压缩的fastq文件 zcat SRR3589956_1.fastq.gz | fastqc -t 4 stdinfastqc SRR3589956_1.fastq.gz结果会得到一个html文件和一个zip压缩包。
绿色表示通过，红色表示未通过，黄色表示不太好。一般而言RNA-Seq数据在sequence deplication levels 未通过是比较正常的。毕竟一个基因会大量表达，会测到很多遍。 具体含义可以看这里： http://jingyan.baidu.com/article/49711c6149e27dfa441b7c34.html 由于有14个结果，如果一个一个打开过去，一定会麻烦死，最好有一种一劳永逸的方法。 利用conda安装软件尤其简单， conda install multiqcmultiqc --help使用也很方便， # 先获取QC结果 ls *gz | while read id; do fastqc -t 4 $id; done# multiqcmultiqc *fastqc.zip --pdf会有一个html文件用来了解总体情况 除了用multiQC查看多个QC结果以外，还可以专门写一个脚本看每个样本的reads数量，GC含量，Q20，Q30的比例 Python脚本 逻辑：
import reimport zipfile# read the zip filedef zipReader(file): qcfile = zipfile.ZipFile(file) data_txt = [file for file in qcfile.namelist() if re.match('.*?_data\.txt', file)][0] data = [bytes.decode(line) for line in qcfile.open(data_txt)] return datadef fastqc_summary(data): module_num = 0 bases = 0 Q20 = 0 Q30 = 0 for line in data: if re.match('Filename', line): filename = line.split(sep='\t')[1].strip() if re.match('Total Sequence', line): read = line.split(sep='\t')[1].strip() if re.match('%GC', line): GC = line.split(sep='\t')[1].strip() if re.match('[^#](.*?\t){6}',line): bases = bases + 1 if float(line.split('\t')[1]) > 30: Q20 = Q20 + 1 Q30 = Q30 + 1 elif float(line.split('\t')[1]) > 20: Q20 = Q20 + 1 if re.match('>>END', line) : module_num = module_num + 1 if module_num >= 2: break Q20 = Q20 / bases Q30 = Q30 / bases summary = [filename, read, GC, str(Q20), str(Q30)] return summaryif __name__ == '__main__': import sys for arg in range(1, len(sys.argv)): data = zipReader(sys.argv[arg]) summary = fastqc_summary(data) with open('summary.txt', 'a') as f: f.write('\t'.join(summary) + '\n') |
|
|
