搜索
查看: 260|回复: 0

[mRNA-seq] 转录组入门(5~9)比对计数

[复制链接]

4

主题

11

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2019-11-17 08:52:20 | 显示全部楼层 |阅读模式
5、数据过滤
尝试删除
6、下载参考基因组下载参考基因组mkdir reference && cd referencemkdir -p index/hg19 && cd index/hg19#wget 下载或者其他工具下载wget http://hgdownload.soe.ucsc.edu/g ... Zips/chromFa.tar.gz tar -zvf chromFa.tar.gzcat *.fa > hg19.farm chr*#axel下载axel http://hgdownload.cse.ucsc.edu/g ... s/chromFa.tar.gztar zvfx chromFa.tar.gzcat *.fa > hg19.farm chr*.facd ~/referencemkdir -p genome/hg38  && cd genome/hg38nohup wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
了解与下载基因组注释gtf文件
基因组标准注释文件-Gencode数据库,由sanger研究所负责整理和维护,主要记录了基因组的功能注释,比如基因组每条染色体上面有哪些编码蛋白的基因,哪些假基因,哪些lncRNA的基因以及它们的坐标,基因上面的外显子内含子坐标,UTR区域坐标。
ftp://ftp.ebi.ac.uk/pub/databases/gencode/
#axel下载axel ftp://ftp.ebi.ac.uk/pub/database ... .annotation.gtf.gz#解压缩gzip -d gencode.v26lift37.annotation.gtf.gz##19年10月份最新版本release 32##axel ftp://ftp.ebi.ac.uk/pub/database ... 7.annotation.gtf.gz
RNA-seq(6): reads计数,合并矩阵并进行注释推荐的小鼠注释文件下载地址:ftp://ftp.sanger.ac.uk/pub/gencode/Gencodemouse/releaseM10/gencode.vM10.chrpatchhapl_scaff.annotation.gtf.gz
#下载最新版本的axel ftp://ftp.ebi.ac.uk/pub/database ... notation.gtf.gzgzip -d gencode.vM23.chr_patch_hapl_scaff.annotation.gtf.gz
7、序列比对下载index文件cd referecemkdir -p index/hisat && cd index/hisat  ##加-p可以直接定义到下级子目录#wget下载 wget -c ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz#axel下载 axel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gzaxel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gzaxel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gztar zxvf hg19.tar.gztar zxvf mm10.tar.gz# 解压得到两个目录 tar -zxvf *.tar.gz # 删除压缩包 rm -rf *.tar.gz
!可忽略!(一般来说用现成的就行了,==把我自己都搞懵了)通过HISAT2的方法进行创建index#电脑配置还行的,或者是没有现成index的,通过HISAT2的方法进行创建#建立基因组+转录组+SNP索引hisat2-build -p 4 genome.fa --snp genome.snp --ss genome.ss --exon genome.exon genome_snp_tranextract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtfhisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19
其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites,但是先做的目的是为了提高运行效率,建立index, 必须选项是基因组所在文件路径和输出的前缀。
fastq格式的reads比对上去得到sam文件
把fastq格式的reads比对上去得到sam文件,接着用samtools把它转为bam文件
一个sam文件16G,心疼我的小电脑……
cd data/191113testmkdir -p aligned hisat2 -t -x "下载的hiast2物种索引文件" -1 *_1.fastq.gz -2 *_2.fastq.gz -S *.sam#全部都使用绝对路径#-x 指定基因组索引#-1 指定第一个fastq文件#-2 指定第二个fastq文件#-S 指定输出的SAM文件###for循环语句##人的样本 humo sapiensfor i in `seq 57 58`dohisat2 -t -x /Users/pipiliaoping/data/reference/index/hisat/hg19/genome -1 /Users/pipiliaoping/data/191113test/fastq/SRR35899${i}_1.fastq.gz -2 /Users/pipiliaoping/data/191113test/fastq/SRR35899${i}_2.fastq.gz -S /Users/pipiliaoping/data/191113test/aligned/SRR35899${i}.samdone##鼠样本 mus musculusfor i in `seq 59 62`dohisat2 -t -x /Users/pipiliaoping/data/reference/index/hisat/mm10/genome -1 /Users/pipiliaoping/data/191113test/fastq/SRR35899${i}_1.fastq.gz -2 /Users/pipiliaoping/data/191113test/fastq/SRR35899${i}_2.fastq.gz -S /Users/pipiliaoping/data/191113test/aligned/SRR35899${i}.samdone
格式转换,排序,索引
Samtools软件,用于处理sam与bam格式的工具软件,能够实现二进制查看、格式转换、排序及合并等功能,结合sam格式中的flag、tag等信息,还可以完成比对结果的统计汇总。同时利用linux中的grep、awk等操作命令,还可以大大扩展samtools的使用范围与功能。
view: BAM-SAM/SAM-BAM 转换和提取部分比对  如sam转bam:samtools view -S ~.sam -b > ~.bam
sort: 比对排序
index: 索引排序比对
SAM格式是目前用来存放大量核酸比对结果信息的通用格式,也是人类能够“直接”阅读的格式类型,而BAM和CRAM是为了方便传输,降低存储压力将SAM进行压缩得到的格式形式。
⚠️版本不同会有不同的表达方式,理解代码及时修改
for i in `seq 59 62`dosamtools view -S SRR35899${i}.sam -b > SRR35899${i}.bamsamtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bamsamtools index SRR35899${i}_sorted.bamdone
8、比对质控(QC)
A survey of best practices for RNA-seq data analysis里面,提到了人类基因组应该有70%~90%的比对率,并且多比对read(multi-mapping reads)数量要少。另外比对在外显子和所比对链(uniformity of read coverage on exons and the mapped strand)的覆盖度要保持一致。
# Python2.7环境下pip install RSeQCbam_stat.py -i SRR3589956_sorted.bamread_distribution.py -i RNA-Seq/aligned/SRR3589956_sorted.bam -r reference/hg19_RefSeq.bed
截取前面100行数据,分别排序然后查看结果
head -1000 SRR3589957.sam > test.samsamtools view -b test.sam > test.bamsamtools view test.bam | head
默认排序是根据染色体位置
samtools sort test.bam defaultsamtools view default.bam | head
Sort alignments by leftmost coordinates, or by read name when -n is used。默认按照染色体位置进行排序,而-n参数则是根据read名进行排序。当然还有一个-t 根据TAG进行排序
samtools sort -n test.bam sort_leftsamtools view sort_left.bam | head
三板斧的view是一个非常实用的子命令,除了之前的格式转换以外,还能进行数据提取和提取。比如说提取1号染色体1234-123456区域的比对read
samtools view SRR3589957_sorted.bam chr1:1234-123456 | head
flagstat看下总体情况
samtools flagstat SRR3589957_sorted.bam
想用samtools筛选恰好配对的read,就需要用0x10
samtools view -b -f 0x10 SRR3589957_sorted.bam chr1:1234-123456  > flag.bamsamtools flagstat flag.bam
9、Reads计数
用htseq对比对产生的bam进行count计数
htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m intersection-nonempty yourfile_name.bam ~/.gtf > counts.txt
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 输出帮助信息。
htseq-count计数得到.count文件
mkdir -p matrix/#for循环的另外一种写法#for i in $(seq 56 58)#htseq-count -r name -f bam ~_sorted.bam 基因组注释gtf文件 > 输出文件~.count#htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m intersection-nonempty yourfile_name.bam ~/reference.gtf > counts.txt
推荐第三种feature count,快到起飞。第二种.log保存程序日志格式
featureCounts比htseq-count快20倍!!!
#人 第一种for i in `seq 56 58`do featureCounts -T 6 -p -t exon -g gene_id -a /Users/pipiliaoping/data/reference/index/gencode.v26lift37.annotation.gtf -o /Users/pipiliaoping/data/191113test/matrix/SRR35899${i}_featureCounts.txt /Users/pipiliaoping/data/191113test/aligned/SRR35899${i}_sorted.bamdone#人 第二种for i in `seq 56 58`dohtseq-count -s no -r name -f bam /Users/pipiliaoping/data/191113test/aligned/SRR35899${i}_sorted.bam /Users/pipiliaoping/data/reference/index/gencode.v26lift37.annotation.gtf 1> /Users/pipiliaoping/data/191113test/matrix/SRR35899${i}.count 2> /Users/pipiliaoping/data/191113test/matrix/SRR35899${i}.logdone#小鼠for i in `seq 59 62`dofeatureCounts -T 6 -p -t exon -g gene_id -a /Users/pipiliaoping/data/reference/index/gencode.vM23.chr_patch_hapl_scaff.annotation.gtf -o /Users/pipiliaoping/data/191113test/matrix/SRR35899${i}_featureCounts.txt /Users/pipiliaoping/data/191113test/aligned/SRR35899${i}_sorted.bam    donecd 路径for i in `seq 56 58`do cut -f 1,7 SRR35899${i}_featureCounts.txt |grep -v '^#' >SRR35899${i}_cut_feacturCounts.txt#基因ID第一行,表达量第七行perl -alne '$sum += $F[1]; END {print $sum}' SRR35899${i}_cut_feacturCounts.txtdone#featureCounts结果总count数perl -alne '$sum += $F[1]; END {print $sum}' feacturCounts.txt#22730409#HTSeq-count结果的总count数grep -v '^_' HTSeq-count.txt |perl -alne '$sum += $F[1]; END {print $sum}'#21514247featureCounts -T 6 -p -t exon -g gene_id -a /Users/pipiliaoping/data/reference/index/gencode.vM23.chr_patch_hapl_scaff.annotation.gtf -o Counts.txt *.bam#这里的*bam文件用大文件的bam和小文件sorted.bam结果一摸一样。multiqc *txt.summary
bam文件出错行(takes from 3 to 5 positional arguments but 6 were given) 导致报错的原因在于htseq不同版本和所处环境
⚠️我的报错:my_showwarning() takes from 3 to 5 positional arguments but 6 were given [Exception type: TypeError, raised in warnings.py:99]
解决办法conda uninstall htseq 卸载htseq旧版本0.7.conda install htseq 重新安装,自动安装最新版本0.11.2
⚠️我的报错:
2617093 GFF lines processed.Warning: Read SRR3589956.268526 claims to have an aligned mate which could not be found in an adjacent line.100000 SAM alignment record pairs processed.Warning: 54028562 reads with missing mate encountered.56101093 SAM alignment pairs processed.
解决办法 (直接换成featureCounts不好吗= =)#htseq-count -r name -f bam ~_sorted.bam 基因组注释gtf文件 > 输出文件~.count-r参数name改为pos即可htseq-count -s no -r pos -f bam ~_sorted.bam 基因组注释gtf文件 > 输出文件~.count
看下每个文件的格式,查看前后10行,第一列ensemblgeneid,第二列read_count计数
$ head -n 10 *.count$ tail -n 10 *.count
QCmultiqc *txt.summary或者multiqc ./






上一篇:转录组入门(1~4)数据下载与QC
下一篇:问问高手,GEO数据处理的论文什么样的期刊接收呀
回复

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-12-6 07:17 , Processed in 0.026054 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.