oppo手机网站建设需求分析,辽宁省建设厅安全员考试官方网站,昆山城乡建设局网站,棋牌网站开发需要多少钱一、序列比对在2016年的一篇综述A survey of best practices for RNA-seq data analysis#xff0c;提到目前有三种RNA数据分析的策略。那个时候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已经被它的作者推荐改用HISAT进行替代。1. Hisat2教程1.1 下载安装#conda直…一、序列比对在2016年的一篇综述A survey of best practices for RNA-seq data analysis提到目前有三种RNA数据分析的策略。那个时候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已经被它的作者推荐改用HISAT进行替代。1. Hisat2教程1.1 下载安装#conda直接安装conda install hisat2#源码下载安装wget wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-source.zipunzip hisat2-2.1.0-source.zipmake1.2 构建index直接下载现有的insex或通过Hisat2的方法进行创建# 其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites但是先做的目的是为了提高运行效率extract_exons.py gencode.v26lift37.annotation.sorted.gtf hg19.exons.gtf extract_splice_sites.py gencode.v26lift37.annotation.gtf hg19.splice_sites.gtf # 建立index 必须选项是基因组所在文件路径和输出的前缀hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg191.3正式比对hisat2基本用法就是hisat2 [options]* -x {-1 -2 | -U } [-S ]基本就是提供index的位置PE数据或者是SE数据存放位置。然而其他可选参数却是进阶的一大名堂。新手就用默认参数呗。hisat2 --dta -p 6 --max-intronlen 5000000 -x Oryza_sativa.IRGSP-1.0.genome -1 C1-1_good_1.fq -2 C1-1_good_2.fq -S C1-1.HISAT_aln.sam hisat2_running.log 211.4 Hisat2输出结果比对之后会输出如下结果解读一下就是全部数据都是100%的2.88%的配对数据一次都没有比对94.20%的数据比是唯一比对2.92%是多个比对。然后如果不按照顺序来有4.96%的比对。之后把剩下的部分用单端数据进行比对的话65.57%数据没比对上33.23%的数据比对一次1.20%比对超过一次。零零总总的加起来是98.20%的比对。20182824 reads; of these:20182824 (100.00%) were paired; of these:581893 (2.88%) aligned concordantly 0 times19011569 (94.20%) aligned concordantly exactly 1 time589362 (2.92%) aligned concordantly 1 times----581893 pairs aligned concordantly 0 times; of these:28886 (4.96%) aligned discordantly 1 time----553007 pairs aligned 0 times concordantly or discordantly; of these:1106014 mates make up the pairs; of these:725197 (65.57%) aligned 0 times367552 (33.23%) aligned exactly 1 time13265 (1.20%) aligned 1 times98.20% overall alignment rate2. SAMtools三板斧SAM(sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式当然他可以用于存放未比对的数据。而目前处理SAM格式的工具主要是SAMToolsview: BAM-SAM/SAM-BAM 转换和提取部分比对sort: 比对排序merge: 聚合多个排序比对index: 索引排序比对faidx: 建立FASTA索引提取部分序列tview: 文本格式查看序列#最常用的三板斧就是格式转换排序索引。samtools view -S SRR35899${i}.sam -b SRR35899${i}.bamsamtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bamsamtools index SRR35899${i}_sorted.bam3. BAM/SAM文件格式SAM文件主要由两个部分构成:header标记了该SAM文件的一些基本信息比如版本、按照什么方式排序的、Reference信息等等。本体每行为一个reads不同列记录了不同的信息列与列之间通过tab分隔。每列的含义MAPQ值表示为mapping的质量值mapping Quality, It equals -10log10Pr{mapping position is wrong}, rounded to the nearest integer, A value 255 indicates that the mapping quality is not available. 该值的计算方法是mapping的错误率的-10log10值之后四舍五入得到的整数如果值为255表示mapping值是不可用的如果是unmapped read则MAPQ为0一般在使用bwa mem或bwa aln(bwa 0.7.12-r1039版本)生成的sam文件第五列为60表示mapping率最高一般结果是这一列的数值是从0到60且0和60这两个数字出现次数最多想把小于2的都丢弃samtools view -bSq 2 file.sam filtered.bamflag的含义1 代表这个序列采用的是PE双端测序2 代表这个序列和参考序列完全匹配没有插入缺失4 代表这个序列没有mapping到参考序列上8 代表这个序列的另一端序列没有比对到参考序列上比如这条序列是R1,它对应的R2端序列没有比对到参考序列上16代表这个序列比对到参考序列的负链上32 代表这个序列对应的另一端序列比对到参考序列的负链上64 代表这个序列是R1端序列 read1;128 : 代表这个序列是R2端序列read2256 代表这个序列不是主要的比对一条序列可能比对到参考序列的多个位置只有一个是首要的比对位置其他都是次要的512 代表这个序列在QC时失败了被过滤不掉了(# 这个标签不常用)1024: 代表这个序列是PCR重复序列(#这个标签不常用)2048: 代表这个序列是补充的比对(#这个标签具体什么意思没搞清楚但是不常用)cigar的含义:cigar中会包含数字代表了特定match持续了多少nt以及不同的字符代表了不同的match情况。30S512M216N12S (30nt soft clip - 512nt exact match - 216nt skipped region - 12nt soft clip)30S (30nt soft clip)40M (40nt exact match)其中不同的字符及其含义如下参考https://www.jianshu.com/p/a584d31418f3https://www.jianshu.com/p/9c87bba244d8二、htseq-count的使用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(质量分析)。具体参考https://www.cnblogs.com/triple-y/p/9338890.htmlhttps://blog.csdn.net/herokoking/article/details/78257714三、基因差异表达分析1. DESeq2(DESeq2不支持无生物学重复的数据)library(DESeq2)#directorydirectorysampleFiles#sampleFilessampleConditionsampleTable#sampleTableddsHTSeq#ddsHTSeqcolData(ddsHTSeq)$conditionddsresres#head(res)png(fileplotMA.png,typecairo)plotMA(dds)dev.off()mcols(res,use.namesTRUE)write.csv(as.data.frame(res),filedeseq2.csv)png(fileplotDispEsts.png,typecairo)plotDispEsts(dds)dev.off()#pheatmapsum(res$padj 0.05, na.rmTRUE)library(pheatmap)select nt log2.norm.counts df pdf(heatmap.pdf,width 6, height 7)pheatmap(log2.norm.counts, cluster_rowsTRUE, show_rownamesFALSE,cluster_colsTRUE, annotation_coldf)dev.off()2. EdgeRcount.matrixID B.count wt.countAT1G01010 384 314AT1G01020 661 861AT1G01030 48 47AT1G01040 5608 6206AT1G01046 77 82AT1G01050 2889 2357AT1G01060 6039 6296AT1G01070 408 240AT1G01073 0 0edgeR.Rlibrary(edgeR)data read.table(count.matrix, headerT, row.names1, com)col_ordering c(1,2)rnaseqMatrix data[,col_ordering]rnaseqMatrix round(rnaseqMatrix)rnaseqMatrix rnaseqMatrix[rowSums(cpm(rnaseqMatrix) 1) 2,]conditions factor(c(rep(B, 1), rep(wt, 1))) #对样本进行分组exp_study DGEList(countsrnaseqMatrix, groupconditions) #建立基因表达列表利用DEGList函数 参数counts为read数矩阵group为上一步的分组变量exp_study calcNormFactors(exp_study) #计算各样本内的标准化因子exp_study_bcv exp_studybcv 0.01 #估计生物学变异系数et exactTest(exp_study_bcv, dispersion bcv ^ 2)#gene decideTestsDGE(et, p.value 0.05, lfc 0)#summary(gene)#head(res)tTags topTags(et,nNULL)result_table tTags$tableresult_table data.frame(sampleAB, sampleBwt, result_table)result_table$logFC -1 * result_table$logFCwrite.csv(result_table,file B-wt.csv)source(Rna-seq/trinityrnaseq-Trinity-v2.4.0/Analysis/DifferentialExpression/R/rnaseq_plot_funcs.R)pdf(B-wt.pdf)plot_MA_and_Volcano(result_table$logCPM, result_table$logFC, result_table$FDR)dev.off()
