查看原文
其他

单细胞测序之原始数据处理

豆豆花花 生信星球 2022-06-07

今天是生信星球陪你的第87天


   你想找辆共享单车,发现满街都是别家车,没有一辆你能骑。

   你想学点生信,搜了“初学者教程”,满眼尽是高大上,没有一句能看懂。

   终于你跨越茫茫宇宙,来到生信星球,发现了初学者的新大陆!


重启文献阅读~

这是单细胞测序第二部分,开始通过真实数据进行数据准备

本次使用的数据是来自2015年Molecular Cell的文献

file:///Users/reedliu1/Library/Mobile%20Documents/com~apple~CloudDocs/223%20-%20Book%20Shelf%20%7C%20%E5%85%A8%E6%96%87%E4%B9%A6%E5%BA%93/Research%20Papers/2015-tech_of_scRNA-Molecular%20Cell.pdf

. 原始数据下载

mkdir raw && cd raw
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR522/ERR522959/ERR522959_1.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR522/ERR522959/ERR522959_2.fastq.gz

. 数据质控

mkdir fastqc && cd fastqc
fastqc -t 10 -o /fastqc ../raw/ERR522959_1.fastq ../raw/ERR522959_2.fastq

结果显示,两个read测试文件用的是Nextera Transposase Sequence的接头,并且都存在污染

原始数据质控

. 数据过滤

使用trim-galore软件,它可以用来切掉存在接头污染或者低质量序列的read末尾,另外过滤后要再用fastqc质控一次

注意:安装时的名称是trim-galore,使用的时候要用trim_galore

mkdir trimmed_fastqc && cd trimmed_fastqc
trim_galore --nextera --fastqc -o ./ ../../raw/ERR522959_1.fastq ../../raw/ERR522959_2.fastq
## --nextera指定污染序列接头类型
## --fastqc在trim完后,设置再次进行fastqc
过滤后的质控

这样的reads才是接下来比对需要的

. 比对须知

比对才能获得测序片段在基因组或转录组上的位点,才能继续分析问题;

要比对reads,就需要有参考基因组和注释文件gtf/gff

  • 对于模式物种,参考基因组和注释文件都可以在三大数据库中下载到:Ensembl、NCBI、UCSC,但不同的是:

  • Ensemble是最容易使用且是最大的注释文件集

  • NCBI只有高质量的基因注释,也就是说,如果研究的物种比较小众,那么NCBI不会入选

  • UCSC也包括大量的基因集、注释集,但标准有点杂

目前做单细胞测序最多的也应该属人和小鼠了,因此参考文件很齐全

  • 简单说明一下GTF文件格式

    (1)seqname:染色体或scaffold编号;

    (2)source:注释来源【比如人类一般采用的都是HAVANA团队注释的】;

    (3)feature:这段区域代表什么基因结构(基因、转录本、外显子?);

    (4)start;(5)end:feature的起始与终止坐标;

    (6)score:feature存在性和它坐标的可信度(可有可无);

    (7)strand:+正链-负链;

    (8)frame:0表示从一个完整地密码子开始,起始于密码子5’端的第一个碱基;1表示在第一个完整密码子前有一个碱基;2表示前面还有2个碱基;而反向链,5’端的第一个碱基在end坐标处;

    (9)attribute:分号分隔的额外信息【可选】包括基因id、转录本id、基因生物类型等等内容

    【内容为空,用.标记】

  • 比对的目的是对基因表达量进行定量或者找出样本间表达差异的基因

  • 使用的软件一:STAR

    它的原理是找出比对到参考基因组中一条或多条序列的最长的测序read
    示意图如下:蓝色的是一条测序read,它横跨了两个外显子;黄色的是可变剪切区域;            
    STAR软件先找到read的第一部分,和参考序列的第一个外显子匹配,另外又发现read第二部分和第二个外显子匹配。它就是用这种方式发现可变剪切事件或者染色体重排的。
    缺点:当参考基因组很大的时候,比如人和小鼠,它是非常占内存的

STAR原理

 第一步:构建索引

mkdir index && cd index    
## vi create_index.sh  
fasta=/DIR to hg19/hg19.fa  
STAR --runThreadN 20 \   
--runMode genomeGenerate \  
--genomeDir ./ \  
--genomeFastaFiles $fasta

【这里用了72线程,共23分钟就跑完了】

第二部:比对

mkdir star_align && cd star_align  
## vi star_align.sh  
STAR --runThreadN 20 \  
--readFilesIn ../trimmed_fastqc/ERR522959_1_trimmed.fq ../trimmed_fastqc/ERR522959_1_trimmed.fq \  
--genomeDir ./ \  
--outFileNamePrefix ./STAR
## output files name prefix (including full or relative path)

  • 使用的软件二:Kallisto

上面用的STAR是将reads直接比对到参考序列,因此比对速度比较慢;而Kallisto是利用k-mer比对参考序列,在转录组测序中也可以用它来分析,它不需要拼接转录本,速度很快,也被称作"pseudo-aligner"

下期预告:什么是k-mer?为何K-mer比直接比对reads速度要快?Kallisto怎么使用?


点击底部的“阅读原文”,获得更好的阅读体验哦😻

初学生信,很荣幸带你迈出第一步。

我们是生信星球,一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言、联系微信或发送邮件到Bioplanet520@outlook.com


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存