Author：ligc
Date：19/5/15
在基因水平上，常用的软件为HTSeq-count，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等。
在转录本水平上，一般常用工具为Cufflinks和它的继任者StringTie， eXpress。
在外显子使用水平上,用于分析差异外显子使用的DEXSeq提供了一个Python脚本（dexseq_prepare_annotation.py）执行这个任务。
使用htseq-count
for i in `seq 59 62`
do
nohup htseq-count -f bam -r pos -s no \
~/AKAP95_rna_seq/alignment/SRR35899${i}.sorted.bam \
~/AKAP95_rna_seq/reference/annotation/gencode.vM21.annotation.gtf \
>SRR35899${i}_matrix.count 2>SRR35899${i}_htseq.log &
done
- -f bam/sam： 指定输入文件格式，默认SAM
- -r name/pos: 你需要利用samtool sort对数据根据read name或者位置进行排序，默认是name
- -s yes/no/reverse: 数据是否来自于strand-specific assay。DNA是双链的，所以需要判断到底来自于哪条链。如果选择了no， 那么每一条read都会跟正义链和反义链进行比较。默认的yes对于双端测序表示第一个read都在同一个链上，第二个read则在另一条链上。
- -a 最低质量， 剔除低于阈值的read
- -m 模式 union（默认）, intersection-strict and intersection-nonempty。一般而言就用默认的，作者也是这样认为的。
-
-i id attribute: 在GTF文件的最后一栏里，会有这个基因的多个命名方式（如下）， RNA-Seq数据分析常用的是gene_id。
htseq_count_result
合并表达矩阵
(AKAP95_workspace) [ligc@cluster reads_count]$ cat merge_read_count.sh
#CMD:python merge_read_count.py EV_3_count.tab,EV_4_count.tab,DNMT3B_2_count.tab,DNMT3B_3_count.tab,DNMT3B_4_count.tab EV_3,EV_4,DNMT3B_2,DNMT3B_3,DNMT3B_4 | less -S
# -*- coding: UTF-8 -*-
import sys
sample_counts = sys.argv[1] ##五个样本
sample_names = sys.argv[2] ##样本名
count_files = sample_counts.split(",")
#print(count_files)
sample_ids = sample_names.split(",")
#print(sample_ids)
count_dict = {} ##字典
for count_file in count_files: ## 循环五个sample，生成gene_id和对应的counts
with open (count_file) as count:
for line in count:
if line.startswith("__"):
continue
line = line.rstrip("\n")
ele = line.split("\t")
#print(ele)
#gene_id -> [count1,count2,count3]
if ele[0] in count_dict:
count_dict[ele[0]].append(ele[1])
else:
count_dict[ele[0]] = [ele[1]]
print("gene_id\t" + "\t".join(sample_ids))
for gene_id in count_dict:
#print(count_dict[gene_id])
print(gene_id + "\t" + "\t".join(count_dict[gene_id]))
merge_count_matrix
R代码
rm(list = ls())
options(stringsAsFactors = FALSE)
raw_count <-read.table(file="merge_count.matrix",sep = "\t",header =T)
library(tidyverse)
new_count_matrix <- separate(raw_count,gene_id,into = "gene_id",sep = "[.]")
summary(new_count_matrix)
GAPDH <- new_count_matrix[new_count_matrix$gene_id =="ENSMUSG00000057666",]
AKAP95 <- new_count_matrix[new_count_matrix$gene_id=="ENSMUSG00000024045",]
