转录组数据分析—HTseq定量

源代码 · 2024-9-2 06:24:22

HTSeq作为一款可以处理处罚高通量数据的python包，由Simon Anders, Paul Theodor Pyl, Wolfgang Huber等人携手推出HTSeq — A Python framework to work with high-throughput sequencing data。自觉布以来就备受广大分析职员青睐，其提供了许多功能给那些认识python的大佬们去自大修改利用，同时也分身着给小白们提供了两个可以拿来可用的可实验文件 htseq-count（计数）和 htseq-qa（质量分析）。

这里必要留意的是HTSeq作为read counts的计数软件，承接的是上游比对软件对于clean data给出的比对结果，而且必须为依据reads名称排序的sort bam文件。
所以在之前需应用samtools对bam文件举行排序。

nohup cat ../SRR_Acc_List.txt | while read line; do samtools sort -o $line\_sorted.bam -n -T sorted $line.bam; done &（依然是批量背景运行）安装

利用conda举行安装 conda install -c bioconda htseq=1.99.2

定量

HTSeq盘算counts数有3种模式，如下图所示（ambiguous表现该read比对到多个gene上；no_feature表现read没有比对到gene上）：

htseq-count对链特异性及非链特异性数据的处理处罚是差别的因此在处理处罚数据前肯定要明确目的数据的建库方式（链特异性或非链特异性的）。明确这一关键题目之后就可以举行定量啦！
#htseq-count常用参数-f | --format    default: sam  设置输入文件的格式，该参数的值可以是sam或bam。-r | --order    default: name  设置sam或bam文件的排序方式，该参数的值可以是name或pos。前者表现按read名举行排序，后者表现按比对的参考基因组位置举行排序。若测序数据是双末了测序，当输入sam/bam文件是按pos方式排序的时间，两端reads的比对结果在sam/bam文件中一样平常不是紧邻的两行，步伐会将reads对的第一个比对结果放入内存，直到读取到另一端read的比对结果。因此，选择pos大概会导致步伐利用较多的内存，它也得当于未排序的sam/bam文件。而pos排序则表现步伐以为双末了测序的reads比对结果在紧邻的两行上，也得当于单端测序的比对结果。许多别的表达量分析软件要求输入的sam/bam文件是按pos排序的，但HTSeq保举利用name排序，且一样平常比对软件的默认输出结果也是按name举行排序的。-s | --stranded    default: yes  设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表现非链特异性测序；假如单端测序，yes表现read比对到了基因的公理链上；假如双末了测序，yes表现read1比对到了基因公理链上，read2比对到基因负义链上；reverse表现双末了测序环境下与yes值相反的结果。根听分析文件的明确，一样平常环境下双末了链特异性测序，该参数的值应该选择reverse（本人临时没有测试该参数）。-a | --a    default: 10  忽略比对质量低于此值的比对结果。在0.5.4版本从前该参数默认值是0。-t | --type    default: exon  步伐会对该指定的feature（gtf/gff文件第三列）举行表达量盘算，而gtf/gff文件中别的的feature都会被忽略。-i | --idattr    default: gene_id  设置feature ID是由gtf/gff文件第9列谁人标签决定的；若gtf/gff文件多行具有雷同的feature ID，则它们来自同一个feature，步伐会盘算这些features的表达量之和赋给相应的feature ID。-m | --mode    default: union  设置表达量盘算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的表现图。从图中可知，对于原核生物，保举利用intersection-strict模式；对于真核生物，保举利用union模式。-o | --samout 输出一个sam文件，该sam文件的比对结果中多了一个XF标签，表现该read比对到了某个feature上。-q | --quiet  不输出步伐运行的状态信息和警告信息。-h | --help  输出资助信息。#我的批量分析代码#写个bashvim htseq.sh cat SRR_Acc_List.txt | while read line do       htseq-count -f bam -s no -m intersection-nonempty 2-mapping/$line\_sorted.bam ../../reference/homo/Homo_sapiens.GRCh38.104.gtf > 3-quantity/$line\_count.txt done#背景运行bashnohup bash htseq.sh &结果如下

接下来我们利用R语言对结果举行归并，先看一下文件内容
##看一下恣意文件前十行head SRR14760842_count.txtENSG00000000003 6157ENSG00000000005 0ENSG00000000419 2925ENSG00000000457 791ENSG00000000460 910ENSG00000000938 0ENSG00000000971 2192ENSG00000001036 3245ENSG00000001084 1427ENSG00000001167 1897#看一下恣意文件后十行tail SRR14760842_count.txtENSG00000288721 6ENSG00000288722 166ENSG00000288723 1ENSG00000288724 0ENSG00000288725 1__no_feature 1980267__ambiguous    1633742__too_low_aQual 821454__not_aligned 961622__alignment_not_unique  1571191观察数据我们可以发现，末了五行是不必要的因此必要去除，此外文件缺少表头信息，所以在处理处罚时要留意标注，同时我们还必要对ensembleid举行解释。上代码！
setwd("F:\\RNA-Seq\\GSE176393（SRP323246）\\3-quantity")library(magrittr)library(clusterProfiler)library(org.Hs.eg.db)library(limma)#读取文件files = list.files(pattern = "_count")file <- list()for (i in files) {  file[]<-read.table(i,sep = "\t",check.names = F,col.names = c("gene_id",i))}rt = Reduce(function(x,y) merge(x,y,by="gene_id",all.x=TRUE),file,accumulate=FALSE)#解释rt <- rt[-1:-5,]ensembl = rt[,1]geneSymbol <- bitr(ensembl, fromType = "ENSEMBL", #fromType为输入数据ID范例                toType = "SYMBOL", #toType是指目的范例                OrgDb = org.Hs.eg.db)#Orgdb是指对应的解释包，这是人的解释包，干系解释包请更换rt <- merge(geneSymbol,rt,by.x="ENSEMBL",by.y="gene_id")#对基因名称重复的reads count举行处理处罚rt<-rt[,-1]rt<-as.matrix(rt)rownames(rt)<-rt[,1]rt<-rt[,2:ncol(rt)]rt<-avereps(rt)rt1<-apply(rt,2,as.numeric)rownames(rt1)<-rownames(rt)rt1<-as.data.frame(rt1)rt1<-cbind(row.names(rt1),rt1)colnames(rt1)[1]<-"genes"#输出write.table(rt1,"rawread.txt",sep = "\t",row.names = F,col.names = T,quote = F)末了就得到了归并的read count文件
genes SRR14760842_count.txt SRR14760843_count.txt SRR14760844_count.txt SRR14760845_count.txt SRR14760846_count.txt SRR14760847_count.txt SRR14760848_count.txt SRR14760849_count.txt SRR14760850_count.txtTSPAN6  6157 5622 6525 4781 4797 4543 3796 3663 4326TNMD 0 0 0 0 0 0 0 0 0DPM1 2925 2901 2703 1977 2176 2214 2166 2153 2416SCYL3 791 764 681 627 575 576 295 336 353C1orf112 910 806 726 527 545 550 521 587 620FGR 0 0 0 1 0 0 0 0 0CFH 2192 2064 2219 1035 958 1036 245 260 317FUCA2 3245 3049 3256 2229 2340 2095 2636 2942 3107GCLC 1427 1428 1495 1390 1509 1281 3693 3543 3891NFYA 1897 1747 1954 1605 1624 1551 1606 1438 1619STPG1 407 395 403 545 540 545 489 421 569NIPAL3  1570 1604 1839 2174 2263 2136 1309 1327 1575LAS1L 2594 2444 2296 1770 1793 1916 1878 2103 2131ENPP4 2015 2023 1998 2121 2189 2065 1867 1625 1916SEMA3F  131 158 136 175 157 165 116 124 193CFTR 310 341 331 161 119 126 39  35  38

转录组数据分析—HTseq定量

所属分类: 问答交流

新帖推荐: 30日

推荐作品