生信技能树学习笔记 比对过程: • 1.建索引 • 2.比对参考基因组 • 3.sam转bam 用到的软件——Hisat2 Hisat2主要是用来进行转录组数据的比对。 sh Hisat2Index.sh >Hisat2Index.sh.log & ## ----比对# 进入比对文件夹cd $HOME/project/Human-16-Asthma-Trans/Mapping 单个样本比对 ## 单个样本比对,步骤分解index=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dnainputdir=$HOME/project 多个样本比对 这里需要用到管道符|串联 比对参考基因组 和 sam转bam两个步骤 这里的2代表下面这个程序中输出的过程,并将其重定向到样本对应的log文件中 关注点: • 总比对率:一般都能在80%以上 *log 结果 可视化结果 比对率过低可能 1.细菌污染 2.核糖体RNA 3.比对文件物种错误 比对结果文件:sam/bam格式 SAM(The Sequence Alignment/Map format
注意,这其中有自身比对的问题。 DNAdot 作者已经不再维护 Mac DNA Strider Mac/Linux,Windows Dotter 和Dotter手册 2 最经典最精确的:动态规划算法 全局比对: 参与比对的两条序列里面的所有字符进行比对 优点是非常精确 缺点是运行时间长,不适合数据量庞大的序列数据库搜索 3 目前大多数数据库搜索工具中使用的算法:BLAST算法(Basic local alignment search tool) 前者适合较少量序列间比对,BLAST适合从一组大量序列中搜索与查询相似的序列 BLAST总体比对算法的思想是:首先通过完全匹配来查找序列,然后通过允许有误匹配的方式来扩展比对区域。 BLAST可以用来做什么 -1 推断和鉴定查询序列的功能 -2 指导实验设计论证该功能 -3 找到在模式生物中与查询序列相似的序列,进一步研究其功能 -4 在目标物种发行与查询序列相似的同原序列
二代高通量测序具有以下特点: 1.测序覆盖全基因组 2.测序数据读长短 3.测序数据具有一定的错误率 4.测序数据深度高 5.测序数据具有 3.2 比对算法 短序列比对有很多比对软件,例如 bwa,soap,bowtie2,hisat2,subread 等,在众多的短序列比对软件中,BWA 几乎已经成为默认的行业标准。 目前 bwa 软件又额外发布了 BWA-MEM2 比对软件,该软件比对速度更快,精确度更高。 bwa-mem2 官网:https://github.com/bwa-mem2/bwa-mem2 3.3 比对结果 pairend 比对 综合考虑两条 reads 与参考序列的比对以及比对错误率情况 1、两条 reads 都比对不上; 2、一条比对上,另外一条比对不上,或者另外一条比对到另外染色体,或者两条比对不在正常 insert size 范围内; 3、一对一比对无错配,
专题目录: 1、第1篇:ATAC-seq的背景介绍以及与ChIP-Seq的异同 2、 这部分内容包括对原始测序数据质控,然后比对过滤,这是所有NGS数据处理的上游分析。 学习目标 用FastQC进行质控检测 用Trimmomatic进行质量过滤 用Bowtie2比对,并理解相关参数含义 测序reads 的质控流程示意图 ? Bowtie2是一个快速精确的比对工具,基于Burrows-Wheeler Transform 构建基因组的FM 索引,比对过程所耗内存少。 Bowtie2支持局部、双端、缺口比对模式,对大于50bp的reads比对效果更好(小于50bp的reads用Bowtie1)。 / Bowtie2 比对 p: 线程数 q: reads是fastq格式 x: index路径 U: fastq路径 S: 输出Sam格式文件 ## 课程中给出的代码是单端比对 bowtie2 -p 2
合并数据前的“风险扫描”:要把两份客户名单合在一起用,得先看看里面有没有同一个客户信息却不一样(比如电话号码不同),不然合并完会更乱。2.谁是“裁判”?比哪些“科目”?这是最关键的一步。 (10000-10200)/10000=-2%,这说明业务系统比财务系统高了2%。看整体统计特征:不一定每条记录都去比。你可以对比两份数据的总和、平均数是不是差不多。 2.问:做模糊匹配时,相似度阈值设成多少合适?90%还是95%?感觉很难把握。答:这个确实没有标准答案,需要结合具体场景来摸索。核心是评估业务风险:你需要思考,相似度低到多少,就会导致错误的业务判断? 3.问:要比对的数据来自两个不同的旧系统,字段名不同、格式混乱、代码值含义也不一样,简直无从下手。答:这是数据比对中最经典、也是最考验耐心的“脏活累活”。 明确系统A里的CUST_CODE字段,对应的就是系统B里的客户编码;系统A用1/2表示性别,系统B用M/F,它们需要互相转换。
·1.参考基因组准备·2.比对:Hisat2 Salmon1.参考基因组准备参考基因组数据库常用参考基因组数据库Ensembl:www.ensembl.org #用得最多数据库完善有基因对应的IDNCBI Hisat2,Subjunc·基因比对:1建索引 2比对参考基因组 3sam转bamHisat2图片----1.构建索引# 进入参考基因组目录cd $HOME/database/GRCh38.105 &----2.比对# 进入比对文件夹cd $HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2## 单个样本比对,步骤分解index=/home/t_rna 2_val_2.fq.gz \ -S ${outdir}/SRR1039510.Hisat_aln.sam #匹配的项目文件比对完生成结果如图a图片----3.sam转bamsamtools )10个样本 转录组估算使用空间:一个样本1.5G大小 *101、质控:cleandata 1.5GG*102、比对: sam 13G10 2(膨胀),bam 2G*10共约 410G简单粗暴 转录组数据多大
2. 运行速度 随着测序价格的下降和数据深入挖掘的需求,测序量越来越大,海量测序reads的比对,要求速度上必须够快。 3. 同时支持DNA和RNA数据的比对,软件官网如下 http://ccb.jhu.edu/software/hisat2/index.shtml 目前最新版为为hisat2. 在进行比对前,首先需要对参考基因组建立索引, 基本用法如下 hisat2-build -p 20 hg19.fa hg19 对于转录组数据,在构建索引时,可以通过gtf文件,得到剪切位点和exon的信息 对于单端数据,采用-U指定输入文件;对于双端数据,采用-1和-2分别指定R1端和R2端的输入文件。 reads比对到基因组上的一个位置,我们称之为一个alignment。 单端数据比对的用法如下 hisat -x hg19 -p 20 -U reads.fq -S align.sam 双端数据用法如下 hisat -x hg19 -p 20 -1 R1.fq -2 R2.
RNA-Seq数据比对和DNA-Seq数据比对有什么差异? RNA-Seq数据分析分为很多种,比如说找差异表达基因或寻找新的可变剪切。 若有多组数据,使用逗号将文件分隔。Reads的长度可以不一致。 -2 <m2> 双端测序结果的第二个文件。若有多组数据,使用逗号将文件分隔,并且文件顺序要和-1参数对应。 -U <r> 单端数据文件。若有多组数据,使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。 多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型,进行比对。 -S <hit> 指定输出的SAM文件。 ,当然他可以用于存放未比对的数据。
从序列数据库搜索,序列拼接到基因蛋白质功能注释,以及进化树构建等,都依赖于分子序列相似性的比较,也就是序列比对。 序列比对的核心作用就是判断是否同源。 blast 比对中默认使用的就是 BLOSUM62 打分矩阵。其中 62 表示用来构建该矩阵的匹配数据集中精确匹配位点要占 62%。 全局比对与局部比对 例如我们现在有两条序列 S1 和 S2,如果采用全局比对,会得到这种比对效果,而采用局部比对,序列中间的 GCG 满足了最优比对。 下载blast数据库 四、blast 数据库 4.1 NCBI blast 数据库 blast 比对需要建立索引,索引 index,是目录的意思。 六、选项参数 blast 常用选项参数 选项 释义 -h 显示选项参数 -help 显示帮助文档 -db 比对数据库 -query 待比对序列 -out 输出文件名 -evalue 比对 e 值 -outfmt
ChIPseq reads 比对在评估读取质量和我们应用的任何读取过滤之后,我们将希望将我们的读取与基因组对齐,以便识别任何基因组位置显示比对读取高于背景的富集。 由于 ChIPseq 读数将与我们的参考基因组连续比对,我们可以使用我们在之前中看到的基因组比对器。生成的 BAM 文件将包含用于进一步分析的对齐序列读取。图片2. 比对4.1. Rsubread我们可以使用 Rsubread 包将 FASTQ 格式的原始序列数据与 mm10 基因组序列的新 FASTA 文件进行比对。 具体来说,我们将使用 align 函数,因为它利用了 subread 基因组比对算法。 .mainChrs"))然后我们可以使用 bowtie2() 函数对齐我们的 FASTQ 数据,指定我们新创建的索引、SAM 输出的所需名称和未压缩的 FASTQ。
生信技能树学习笔记 subread 官网:http://subread.sourceforge.net/ 构建索引: subjunc:subread-buildindex 5款流行比对工具大比拼:https | while read iddo subjunc -T 5 -i ${index} -r ${inputdir}/${id}_1_val_1.fq.gz -R ${inputdir}/${id}_2_ val_2.fq.gz -o ${outdir}/${id}.Subjunc.bam 1>${outdir}/${id}.Subjunc.log 2>&1 && samtools sort -@ 6 - Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam ##----depth统计测序深度# 得到的结果中,一共有3列以指标分隔符分隔的数据 reads,flag值的理解# (0x100) 代表着多比对情况,所以直接用samtools view -f 0x100可以提取 multiple比对的 情况
localhost:19200/index01/tweet' 'http://localhost:19200/index01/tweet'Unchanged 1Unchanged 2Unchanged http://localhost:19200/index01/tweet' 'http://localhost:29200/index01/_doc'Unchanged 1Deleted 2Updated http://localhost:19200/index01/tweet' 'http://localhost:39200/index01/_doc'Unchanged 1Deleted 2Updated
STAR是一款RNA_seq数据专用的比对软件,比对速度非常快,最大的优势是灵敏度高,GATK推荐采用STAR比对,然后进行下游的SNP分析。 单端数据比对的基本用法如下 STAR \ --runThreadN 20 \ --genomeDir hg19_STAR_db \ --readFilesIn reads.fq \ --sjdbGTFfile hg19.gtf \ --sjdbOverhang 149 \ --outFileNamePrefix sampleA \ --outSAMtype BAM SortedByCoordinate 双端数据比对的基本用法如下 ,STAR官方更推荐使用2-pass比对模式,即比对两次,有以下两种方式 multi-sample 2-pass 第一次比对和上述的用法一致,比对完之后,每个样本会产生一个intron的区间文件SJ.out.tab per-sample 2-pass 对于单个样本,在比对时直接添加--twopassMode Basic参数,软件会自动进行两次比对,将第一次比对的SJ.out.tab加入到索引,然后重新比对。
上一篇文章双序列比对与BLAST介绍了两条序列之间进行比对的算法原理及其实现方法,双序列比对常用于同源分析、蛋白质结构推断、相似片段搜寻与数据库比对检索、基因注释等。 考虑m个长度为n的序列进行多序列比对,首先基于二维的推广其算法复杂度有一个O(nm)项;其次在关于最优得分S的推导时其有m个序列则有2m-1个前导项,例如对于三序列(u、v、w)比对有以下7个前导项: 因此标准动态规划法多序列比对的算法复杂度高达O(nm2m),这使得标准动态规划算法进行大量较长的多序列比对十分困难。 根据基准测试数据的研究基于一致性方法的多序列比对产生的结果经常比渐进多序列比对更加准确。 该软件参数众多,但提供了精确度不同的三个常用模式,以适用不同数据集大小、序列保守性的场景: mafft --maxiterate 1000 --localpair in > out #最准确的方法,
全局比对与局部比对有什么不同呢。全局序列比对尝试找到两个完整的序列之间的最佳比对。而局部序列比对不必对两个完整的序列进行比对;可以在每个序列中使用某些部分来获得最大得分。 例如我们现在有两条序列 S1 和 S2,如果采用全局比对,会得到这种比对效果,而采用局部比对,序列中间的 GCG 满足了最优比对。 ,对资源的消耗比较少,官方的给出的数据是两个 5M 左右的基因组,只用 20 秒左右的时间就可以比对完成,消耗的内存大约是 90M,它是使用一种后缀树的算法。 Mummer 官网介绍该软件是一个多才多艺的软件包,因为它可以完成生物数据分析中很多的功能。Mummer 其实是一个软件包,里面包含了很多工具,这些工具搭配起来使用,可以完成非常多的工作。 .filter -l 10000 >kmer45.tiling #mummerplot绘图 mummerplot -p p1 nucmer.filter --png mummerplot -p p2
ChIPseq reads 比对 在评估读取质量和我们应用的任何读取过滤之后,我们将希望将我们的读取与基因组对齐,以便识别任何基因组位置显示比对读取高于背景的富集。 由于 ChIPseq 读数将与我们的参考基因组连续比对,我们可以使用我们在之前中看到的基因组比对器。生成的 BAM 文件将包含用于进一步分析的对齐序列读取。 2. 比对 4.1. Rsubread 我们可以使用 Rsubread 包将 FASTQ 格式的原始序列数据与 mm10 基因组序列的新 FASTA 文件进行比对。 具体来说,我们将使用 align 函数,因为它利用了 subread 基因组比对算法。 .mainChrs")) 然后我们可以使用 bowtie2() 函数对齐我们的 FASTQ 数据,指定我们新创建的索引、SAM 输出的所需名称和未压缩的 FASTQ。
如图,左侧 10 条数据是先前下载的,右侧少了 1 条(数据是随便编的): ? ? Python 操作 因为对 Excel 的函数操作不太熟,第一时间我是用 Python 来比对数据的:选取两份表格中的 id 列,分别复制到两份 txt 文档中,转化为 Python 读取 txt 文档数据 f.readlines() # data2 为 ['9\n', '127\n', '73\n', '44\n', '20\n', '96\n', '1\n', '28\n', '12'] # 对读取到的数据做下简单处理 ,去掉字符串中的换行符 data1 = [x.strip() for x in data1 ] data2 = [y.strip() for y in data2 ] # 选取在 data1 中出现过 、但 data2 中却不包含的数据 result = [i for i in data1 if i not in data2 ] print(result) # 得到结果 ['5'] 根据得到的结果 5
我前面已经对数据进行了质控: 转录组分析 | fastqc进行质控与结果解读 转录组分析 | 使用trim-galore去除低质量的reads和adaptor 接下来我们进行序列比对,利用的软件是 2.采用了新的比对策略 RNA-seq产生的reads可能跨长度比较大的内含子,哺乳动物中甚至最长能达到1MB,同时外显子比较短,read也比较短,会有很多read(模拟数据中大概34%)跨两个外显子的情况 使用Hisat2进行序列比对 创建输出数据的文件夹 mkdir cleandata/hisat2_mm10data 因为比对这一过程很耗内存,所以样本多话,计算机内存不够大,需要分批比对,我就先介绍一个样本比对的命令 /cleandata/hisat2_mm10data/ 原本22G的数据,运行结束后155G,所以自己用本地计算机跑的话,数据多注意一下内存。 ? 四.结果解读 我们比对后得到的是sam文件,关于sam文件是一个什么样的文件,参考文章:生信中常见的数据文件格式。
我前面已经对数据进行了质控: 转录组分析 | fastqc进行质控与结果解读 转录组分析 | 使用trim-galore去除低质量的reads和adaptor 接下来我们进行序列比对,利用的软件是Hisat2 2.采用了新的比对策略RNA-seq产生的reads可能跨长度比较大的内含子,哺乳动物中甚至最长能达到1MB,同时外显子比较短,read也比较短,会有很多read(模拟数据中大概34%)跨两个外显子的情况 使用Hisat2进行序列比对创建输出数据的文件夹 代码语言:javascript代码运行次数:0运行AI代码解释mkdir cleandata/hisat2_mm10data因为比对这一过程很耗内存,所以样本多话 /cleandata/hisat2_mm10data/原本22G的数据,运行结束后155G,所以自己用本地计算机跑的话,数据多注意一下内存。 四.结果解读 我们比对后得到的是sam文件,关于sam文件是一个什么样的文件,参考文章:生信中常见的数据文件格式。
这个专题分享点日常运维中用到的Python脚本 在做数据库迁移后,我们可能需要知道我们的表,索引,存储过程等对象是否迁移成功 这时可以用如下脚本来进行检查 ---- 环境准备 操作系统: Windows ---- 2. 保存源库和目标库信息至文件 将源端的的信息保存在migrate_from.txt文件中 将目标库的信息保存在migrate_to.txt文件中 注意:最后一行时不要换行 ? 脚本内容 #coding=utf8 import os diff1=[] diff2=[] migrate_from=open(r"c:\migrate_from.txt") migrate_to=open for i in migrate_to: diff1.append(i.strip().replace('\t',' ')) for k in migrate_from: diff2. print(j+'\n') print ("迁移后多余的对象\n") for l in diff1: if l not in diff2: print( l+'\n')