生信马拉松 Day16 转录组RNA-seq-1
原创
今天开始学习转录组测序的上下游分析了，艰苦地学完Linux之后感觉生活又光明了！虽然Linux的命令还不太熟悉，只能先想要做什么再从笔记里搜搜关键词，但基本上能跟得上，不用像学Linux一样漏一分钟就要反复回放好多遍才听懂
快乐[]~(￣▽￣)~*
转录组学的用途
1、识别特定条件下上调或下调的基因
2、理解疾病的状态下哪些基因发生改变
3、发现新的基因和转录本：新的可变剪切
转录组学是在整体水平上研究细胞中基因转录情况及调控规律的学科，从RNA水平研究基因的表达情况
对象：某一种特定细胞及组织
测定的结果：几乎所有的转录本及基因序列，狭义上指信使RNA（mRNA）
依赖二代测序
主要用来研究基因表达量、基因功能
信使RNA含量非常少，大概只有5%
上游：Raw data——质控——比对——定量，基于Linux基础
下游分析：差异表达分析——差异基因功能注释+差异基因功能富集，基于R语言
上机测序完成后得到的测序数据：FASTQ文件
一般以fastq或fq结尾保存在FASTQ格式的文件中
下载工作：
SRA/ENA/GEO在线工具
curl/wget/aria2/axel：通用下载方式
prefetch：NCBI提供的SRA数据下载工具
pysradb (Python)：SRA/ENA/GEO多功能下载助手
fastq
高通量测序，荧光发射得到的原始图像数据文件，因此可能有0，经碱基识别分析转化为原始测序序列，称为Raw Data或Raw Reads，结果以FATQ文件格式存储，包含测序序列信息ATCG以及测序质量（指质量好坏，可信度）
文件每个Read由四行描述
1.@开头，之后为序列的标识符及描述信息
2.序列信息，如ATCG，和fa不同，不能写不完就折到下一行
3.以+开头，之后可以再次加上序列的标识及描述信息（保留行），新的数据一般只有+，老数据可能有其他字符
4.碱基质量打分的ASCII码，与第二行的序列相对应，长度必须相同，否则是损坏的文件
碱基质量值（Quality Score或Q-score）是碱基识别（Base Calling）出错的概率的整数映射
例如质量值为10，出错概率就是10%，识别精度就是90%，碱基质量是20，识别出错概率就是1%，识别精度就是99%，碱基质量值为0.10%，那么碱基识别精度就是99.90%，常见标识为Q20，Q30，因此碱基质量值越高越好，说明准确度高
质量的取值一般是0-41
ASCII码前面的数字不方便使用，所以会采用+33或+64来处理
+33取值空间是！~J，+64的取值范围为@~j
因此看到！等符号，就是采用+33的，看到小写字母，就是+64
不同测序仪和测序平台采用的方法不同
到illumina 1.8之后测序仪都是+33，古早的测序数据可能会有+64
对示例文件SRR1039510_1.fastq.gz进行处理
#1.fastq文件统计reads，6种方法
zless SRR1039510_1.fastq.gz | grep "@SRR" |wc -l
zless SRR1039510_1.fastq.gz | grep '^@SRR' |wc -l
zless -S SRR1039510_1.fastq.gz | paste - - - - |wc -l #这个方法的paste - - - -是把文件转换为4列，此时统计有多少行就可以知道有多少reads
zless SRR1039510_1.fastq.gz |wc -l | awk '{print $0/4}'
zless -S SRR1039510_1.fastq.gz |awk '{ if(NR%4==2) {print} }' |wc -l
zless SRR1039510_1.fastq.gz |sed -n '1~4p' |wc -l
#2.输出fastq文件的序列ID（即第一行），4种方式
zless SRR1039510_1.fastq.gz | grep '^@SRR' |less -S
zless SRR1039510_1.fastq.gz | paste - - - - |cut -f 1 |less -S
zless -S SRR1039510_1.fastq.gz |awk '{if(NR%4==1){print}}' |less -S
less fastq\_25000/SRR1039510\_1.fastq.gz | sed -n '1~4p' | less -SN
#3.输出fastq文件的序列ID（即第二行），3种方式
zless SRR1039510_1.fastq.gz | paste - - - - |cut -f 2 |less -S
zless SRR1039510_1.fastq.gz |awk '{if(NR%4==2){print}}' |less -S
less fastq\_25000/SRR1039510\_1.fastq.gz | sed -n '2~4p' | less -SN
#4.统计碱基总数，3种方式
zless SRR1039510_1.fastq.gz |paste - - - - |cut -f 2 |tr -d '\n' |wc -m
zless -S SRR1039510_1.fastq.gz |paste - - - - |cut -f 2 |grep -o [ATCGN] |wc -l
zless -S SRR1039510_1.fastq.gz |awk '{ if(NR%4==2){print} }' | awk 'BEGIN {num=0} {num=num+length($0)} END{ print "num="num}'
#5.查看reads拷贝数
zless SRR1039510_1.fastq.gz | paste - - - - | cut -f2 | sort | uniq -c | sort -nr | head
自己试着玩了一下代码
zless SRR1039510_1.fastq.gz | paste - - - - | cut -f2 | sort | uniq -c | sort -nr | wc -l
#24482
#看看不重复的拷贝数有多少个
zless SRR1039510_1.fastq.gz | paste - - - - | cut -f2 | sort | uniq -c | sort -nr | awk '{print $1}' | awk 'BEGIN {num=0} {num=num+$0} END{ print "num="num}'
#25000
#以极为复杂的方式完成了前面1里的reads值统计原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。