生物信息
记录及分享生物信息相关软件、算法等
Q.1
啥都不会
展开
-
C++下的序列比对NeedlemanWunsch
最近使用C++写了NeedlemanWunsch算法,算法说明请参考之前用python实现的算法文章《基于numpy的基因序列比对算法NeedlemanWunsch》这里就直接贴代码了:#include <iostream>using namespace std;string target_str="ATCGATCGTG";string query_str="ACGACTACGTG";const int match_score = 1;const int mismatch_s原创 2021-09-03 16:10:57 · 637 阅读 · 0 评论 -
用python完成二代测序捕获区间设计
panel的设计其实很简单,根据实验目的来选择需要捕获的区域,我们需要做的就是把这些需要捕获的区域做成一个bed文件。下面就以BRCA1/2两个基因来举例子,一般bed都是设计在基因的CDS区,因为内含子区域往往包含很多低复杂度区域(比如重复区域),所以内含子的捕获性能往往较差,后期分析难度也高。我们需要先准备基因组注释文件,我从NCBI下载的最新版gtf文件(https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh37_latest/原创 2020-12-24 17:01:33 · 667 阅读 · 4 评论 -
病原菌基因组快速搜索算法实现
算法来自(Wellcome Trust Centre for Human Genetics, University of Oxford)19年发表在NBT上的一篇文章这是文章的算法的示意图以及文章中的说明这个算法主要基于布隆过滤器(BloomFilter),如文章中所示,首先我们要拟定几种不同的hash生成器,但是需要注意不能用python的hash()模块,因为hash()在每次重新调用脚本时生成的值都不一样。import hashlibdef my_hash(k_mer, index_原创 2020-06-23 15:01:35 · 367 阅读 · 0 评论 -
将cram/bam文件转换为fastq文件
NCBI下载的cram文件无法直接使用,需要先转成bam/sam文件,根据官网说明下载了cramtools,发现早已没有维护,报错如下:$ java -jar cramtools-3.0.jarError: Invalid or corrupt jarfile cramtools-3.0.jar所以就直接用samtools来转换,但是直接转换会报错:$ samtools view -b NA12878.final.cram > NA12878.bam &Failed to popu原创 2020-06-19 15:58:14 · 3461 阅读 · 0 评论 -
使用gatk检测WES数据中的cnv
gatk的cnv流程对环境依赖较高,需要调用许多python包,推荐在dockerhub里找官方镜像,或者用conda来配置环境。1、dockerhub 在本地的docker环境中直接拉取镜像,如果没有root权限就用conda安装。docker pull broadinstitute/gatk:4.1.6.02、conda 先下载一个miniconda或者anaconda,然后下载好ga...原创 2020-04-26 16:18:26 · 6494 阅读 · 12 评论 -
关于生信常用的PCA算法
PCA是一种常用的数据降维算法,主要就是根据特征值提出特征值较大的几个特征向量,并将高维向量投影到特征向量上以达到数据降维的目的。第一步当然是处理数据,将样本数据按列放入numpy的数据框(比如100个样本,每个样本4个数据,就是4行100列,反过来也可以,不过要注意调换内积的顺序)。数据需要先做归一化,以减去所有样本的均值来实现。from sklearn.datasets import l...原创 2020-01-09 18:06:12 · 1875 阅读 · 0 评论 -
基于numpy的基因序列比对算法NeedlemanWunsch
这是基于动态规划的一种文本比对算法,常用于基因组序列的比对。比如基因组上有一段序列为ATCGATCTGT,我需要比对的序列为ATCCATCAG,那么首先需要构建一个初始的打分矩阵: A. T. C. G. A. T. C. T. G. T. [[ 0. -1. -2. -3. -4. -5. -6....原创 2020-01-09 16:16:42 · 1138 阅读 · 1 评论 -
通过贝叶斯公式计算阳性概率
今天有一个实习的小朋友说他们课程有一道题不知道怎么做,需要用贝叶斯公式计算阳性概率,题目如下:已知某病发病率5%,假阴率5%,假阳率5%,如果一个人检出了阳性结果,那么他患该病的可能性有多大?首先需要说明,假阳就是把没病的检成有病的,假阴就是把有病的检成没病的。然后搬出贝叶斯公式:p(x∣y)=p(x)∗p(y∣x)p(y)p(x|y) = \frac {p(x)*p(y|x)} {...原创 2020-01-07 16:16:08 · 7505 阅读 · 0 评论 -
python中的T检验
在python中,检验常常是调用scipy.stats进行计算(官方文档:https://docs.scipy.org/doc/scipy/reference/stats.html#module-scipy.stats)T检验也是分为三大块:单样本T检验独立样本T检验配对样本T检验1、先说单样本T检验,我们先通过numpy生成一些正态分布的数据,然后看下单样本T检验from sci...原创 2020-01-03 16:16:44 · 10891 阅读 · 0 评论 -
使用bcftools合并vcf过程中遇到的问题与解决方案
首先因为我的vcf并不是标准的vcf,我用的这些vcf是通过python拼凑了几个call snp软件的结果,所以遇到了各种问题,开帖记录一下。要对多个vcf文件进行合并,就是标准的压缩,做索引,然后merge:bcftools view sample.vcf -Oz -o sample.vcf.gzbcftools index sample.vcf.gzbcftools merge sa...原创 2019-12-27 15:39:22 · 9994 阅读 · 6 评论 -
sra文件下载及解析的问题
昨天在NCBI上下载了一份sra文件,使用fastq-dump提取fastq文件时遇到报错:2019-12-24T08:59:08 fastq-dump.2.9.6 sys: timeout exhausted while reading file within network system module - mbedtls_ssl_read returned -76 ( NET - Readi...原创 2019-12-24 17:21:27 · 4480 阅读 · 0 评论