搜索
查看: 3245|回复: 2

[mRNA-seq] 转录组入门分析6_reads计数

[复制链接]

9

主题

15

帖子

130

积分

注册会员

Rank: 2

积分
130
发表于 2017-8-30 10:44:37 | 显示全部楼层 |阅读模式
本帖最后由 laofuzi 于 2017-8-30 10:44 编辑

1. 要求
1.1用htseq-count分析表达量。
1.2需要用脚本合并所有的样本为表达矩阵。
1.3 对这个表达矩阵可以自己简单在excel或者R里面摸索,求平均值,方差。
1.4 一些生物学意义特殊的基因表现如何,比如GAPDH,β-ACTIN等等。

计数分为三个水平: gene-level, transcript-level, exon-usage-level。在基因水平上,常用的软件为HTSeq-count,featureCounts,BEDTools,Qualimap, Rsubread, GenomicRanges等。在转录本水平上,一般常用工具为Cufflinks和它的继任者StringTie, eXpress。这些软件要处理的难题就时转录本亚型(isoforms)之间通常是有重叠的,当二代测序读长低于转录本长度时,如何进行区分?这些工具大多采用的都是expectationmaximization(EM)。在外显子使用水平上,其实和基因水平的统计类似。但是值得注意的是为了更好的计数,我们需要提供无重叠的外显子区域的gtf文件[2]。用于分析差异外显子使用的DEXSeq提供了一个Python脚本(dexseq_prepare_annotation.py)执行这个任务。
标准化方法: FPKM RPKM TMM TPM。但是,有一些下游分析的软件会要求是输入的count matrix是原始数据,未经标准化,比如说DESeq2,这个时候你需要注意你上一步所用软件会不会进行标准化。
2.重新对bam文件进行排序
#在对双端数据进行处理时,必须要按reads name进行排序,前期是按照默认的染色体位置进行排序的,所以要重新进行排序。编辑如下脚本,另存为sam_name.sh
[AppleScript] 纯文本查看 复制代码
[/color][/align][align=left][color=#000000]#!/bin/bash
[/color][/align][align=left][color=#000000]#转到特定文件夹
cd /mnt/d/AKAP95/aligned
# 将所有的bam文件进行排序
# 不能对reads name 排序的bam进行index
for i in `seq 56 62`
  do samtools -n sort SRR35899${i}.bam -o SRR35899${i}.nsorted.bam
done

sh /mnt/d/AKAP95/aligned/sam_name.sh


3. 下载注释文件
[AppleScript] 纯文本查看 复制代码
[/align][align=left]#注释文件的准备
cd /mnt/d/reference/genomes/hg19 
#下载gtf
wget [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_26/GRCh37_mapping/gencode.v26lift37.annotation.gtf.gz]ftp://ftp.sanger.ac.uk/pub/genco ... 7.annotation.gtf.gz[/url]
#解压Gtf文件, -d 将压缩文件解压。
gzip -d gencode.v26lift37.annotation.gtf.gz

mkdir /mnt/d/reference/genomes/mm10
cd /mnt/d/reference/genomes/mm10
axel [url=ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_M10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz]ftp://ftp.sanger.ac.uk/pub/genco ... f.annotation.gtf.gz[/url]
gzip -d gencode.vM10.chr_patch_hapl_scaff.annotation.gtf.gz


4.htseq-count计数
[Shell] 纯文本查看 复制代码
[/align][align=left]#创建文件夹
mkdir /mnt/d/AKAP95/matrix 

#!/bin/bash
cd /mnt/d/AKAP95/matrix
for i in `seq 56 58`
do
htseq-count -s no -f bam /mnt/d/AKAP95/aligned/SRR35899${i}.nsorted.bam  /mnt/d/reference/genomes/hg19/gencode.v26lift37.annotation.gtf  >/mnt/d/AKAP95/matrix/SRR35899${i}.count  2>/mnt/d/AKAP95/matrix/SRR35899${i}.log
done

for i in `seq 59 62`
do
htseq-count -s no -f bam /mnt/d/AKAP95/aligned/SRR35899${i}.nsorted.bam  /mnt/d/reference/genomes/mm10/gencode.vM10.chr_patch_hapl_scaff.annotation.gtf  >/mnt/d/AKAP95/matrix/SRR35899${i}.count  2> /mnt/d/AKAP95/matrix/SRR35899${i}.log
done

#将脚本存于matrix目录下,命名为read_count.sh。
#运行该脚本
sh /mnt/d/AKAP95/matrix/read_count.sh


Tips:count文件的格式:基因ID一列+reads计数一列

Tips: htseq-count 可以接受pos排序的文件和name排序的文件,按read name重新进行排序,后续reads count会更快,出错率低。


5.构建表达矩阵
[AppleScript] 纯文本查看 复制代码
[/align][align=left]#使用R将这四个文件进行合并,得到最后的融合表达矩阵。
#导入数据,合并表达矩阵
# 分别读取7个文件
options(stringsAsFactors = FALSE)
h_control <- read.table("/mnt/d/AKAP95/matrix/SRR3589956.count", sep="\t", col.names = c("gene_id", "h_control_293"))
h_case_1 <- read.table("/mnt/d/AKAP95/matrix/SRR3589957.count", sep="\t", col.names = c("gene_id","h_miR8"))
h_case_2 <- read.table("/mnt/d/AKAP95/matrix/SRR3589958.count", sep="\t", col.names = c("gene_id","h_miR12"))
m_control_1 <- read.table("/mnt/d/AKAP95/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","m_control1"))
m_control_2 <- read.table("/mnt/d/AKAP95/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","m_control2")) 
m_ case_1  <- read.table("/mnt/d/AKAP95/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap95_1")) 
m_ case_2  <- read.table("/mnt/d/AKAP95/matrix/SRR3589962.count", sep="\t", col.names = c("gene_id","akap95_2"))

#数据整合
# 将三个人类数据和四个小鼠的矩阵分别按照gene_id进行合并,并赋值给各自的raw_count
h_raw_count <- merge(merge(h_control, h_case_1, by="gene_id"), h_case_2, by="gene_id")
m_raw_count <- merge(merge(m_control_1, m_control_2, by="gene_id"), merge(m_case_1, m_case_2, by="gene_id"))
write.csv(h_raw_count, file = "/mnt/d/AKAP95/matrix/h_raw_count.csv", quote = FALSE, row.names = FALSE)
write.csv(m_raw_count, file = "/mnt/d/AKAP95/matrix/m_raw_count.csv", quote = FALSE, row.names = FALSE)
#对表达矩阵进行简单分析
summary(h_raw_count)
summary(m_raw_count)
6. 查看GAPDH、beta-actin基因的Reads count

#human GAPDH
#human beta-actin
#mouse GAPDH
# mouse beta-actin

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:为什么我无法从lncRNA注释文件中获得lncRNA转录本
下一篇:耳聋基因检测数据
回复

使用道具 举报

3

主题

20

帖子

156

积分

注册会员

Rank: 2

积分
156
发表于 2017-10-10 21:41:11 | 显示全部楼层
厉害了.     
回复 支持 反对

使用道具 举报

0

主题

2

帖子

73

积分

注册会员

Rank: 2

积分
73
发表于 2018-8-4 12:01:51 | 显示全部楼层
写的挺仔细的,过来学习学习。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-5-21 15:47 , Processed in 0.042746 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.