浮生终有醒

这个博客没有博主了这个博客没有博主了这个博客没有博主了。。。

1. codingGene的结构codingGene一般表现出来至少是四个结构：TU，MODEL，EXON，CDS。TU不必说了，gene的转录单元；MODEL比较有意思，其实就是转录本的不同剪切方式；EXON是包含UTR部分，CDS不包含UTR。对于noncodingGene，没有CDS的说法，也没有UTR的说法，就不必细说了。 2. 3`UTR的结构3`UTR是PolyA尾

Duplicate是个老大难问题，但处理与否要看具体情况，比如做DNA样本的时候，一定会处理，而RNA样本选择不处理。接下来有一些解决方法，但是“但是”也会很多，接受现实吧~~~首先Duplicate出现的类型有两种，一种是由于PCR扩增的原因导致的完全一样的reads，另一种是比对到基因组上同一位置不同的reads，但由于质量问题、测序错误、比对错误、等位基因等等，被认为是Duplicate

之前看过前辈用python转换矩阵，但python一直没系统学过，所以从网络中学到perl的矩阵转换，两个做了一下比较：import sysfile = open(sys.argv[1], 'r')arr = []for line in file: info = line.rstrip().split() arr.append(info) tarr = [[r

二项分布是重复n次的实验，且每次实验都是独立的，只有两种结果，并且相互对立的，生活中最常见的是投硬币~~~在生物领域内也有很多符合此类分布的，如二倍体动物等位基因，来源于父本和母本的重组等。具体公式什么的博主就不写了，写个关于ASE的例子吧。# cat binom.r | R --slave --args args <- commandArgs()fa <- read.table(args

最开始画频数分布图用的是excel的数据透视表，虽然简单易行，但是画的图杂志不认可，没办法就到处请教别人R的画法。fa=read.table("fa.rrr",header=FALSE)ea=read.table("ea.rrr",header=FALSE)mat_fa=density(data.matrix(fa))#频数分布mat_ea=density(data.matrix(ea)

博主在做PCA的时候，经常遇到2D图无法区分开各个条件的样本，而3D却能直观地感受到样本之间的距离与差异。使用软件是EPD中的IDLE，输入文件为3维的PCAscore结果，参考之前的PCA脚本即可，3D脚本如下：#!pythonimport sys, reimport numpy as npimport matplotlib.pyplot as pltfrom mpl_toolkit

哥们找个图让我画，查了一番和妹子的点拨，总结了一些东西：fa<-read.table("KEYNAME.cluNum.ladder",header=FALSE)mat <- fa[,2]pdf("bar.pdf")barplot(mat,ylim=c(0,7000),width=1,space=0)axis(1,1:96,labels=fa$V1,tick=FALSE)color <

现在根据上一篇对染色体的控制，这篇讲的是染色体自己内部的控制，主配置文件和breaks文件如下：>>>>karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = nochromosomes = hs

karyotype核型图，一般的核型图是直线模式，circos将其转换成圆形。karyotype文件格式分两部分，第一部分为染色体的总长和标识，第二部分为每条染色体的基因区域，当然第一部分是必要的。看了一下这一节的内容发现跟第一节的内容有些不一样，看来必须重新写完整的了。#ideogram.confdefault = 0.0025r#break = 0.5r 同一染色体内部分开

怎么还是helloword？？？！！！事实上，这一章讲的是ideogram的一些参数问题，大体的讲了一些参数所覆盖的情况，所以只是helloworld。本章和第一章内容差不多，先将整体基础配置文件给出，然后一节一节地添加参数画图，最后加入自己的测试文档，搞定~~~#bands.confshow_bands = yesfill_bands = ye

跟上篇几乎一样，主配置文件有变化，其他无变化：karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = nochromosomes = /hs[1234]$/chr1* = redchr2* =

继续翻译，这一节是染色体的分割模式，两个次要配置文件没有变化（ticks和ideograms），主文件变化为：karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = no #不让其显示默认的24条染色体chromosomes

karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000 #默认染色体长度单位为MB>> #类似java类函数的调用> > >以上为主文件，现在看看调用里面是什么吧：show_ticks = yes #刻度尺show_ti

根据网上其他人的建议，符合博主生物类的风格，其实胡说的- -！# .bashrc# Source global definitionsif [ -f /etc/bashrc ]; then . /etc/bashrcfi# User specific aliases and functionsalias le="less -SN";alias l="ls -lhrt";al

在刚开始学的时候不会用R来计算相关性系数，也不会画图，结果博主很悲催地用perl的svg进行画图，很久之前的作品：#!/usr/bin/env perluse warnings;use strict;use lib "/bin/svg_lib";use PLOT qw(Paper End Point Rect Line Text Polyline Path);die "Usage:

生物信息里面也有进制的转换关系，不多说，把简单的写一下：sub d2b()#2to10{ my $bin = shift; my $mod = $bin % 2; return $bin if $bin < 2; $bin = ($bin - $mod) / 2; return &d2b($bin).$mod;}sub dec2bin()#2to10{ my $dec =

heatmap其实博主对这章有一些细节还没有完全理解，照着画倒是可以完成，但是就是心里挺不安的。。。karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = nochromosomes = /hs[1234]

axe和background属于坐标轴，前者是细线条，后者是宽线条- -！！！能这样翻译么- -！！！按照字面意思来讲，background是背景的意思，说那么多还不如画出来实在，例子添加了两个conf，但为了大家能看明白，将其写在主配置文件内比较容易理解，当然也可以多写个conf，使主配置文件简洁一些：karyotype = data/karyotype/karyotype.human.tx

练习一下大名鼎鼎的circos，画图风头正劲，各种CNNS引用。话说是子江兄引进国内，但是从安装到调试，从语法到逻辑，真心说比较复杂，不过看到结果的话心理面就好受点。。。这个博主不知道是算原创还是算翻译，不过由于大部分都是自己在做测试，顺便就算翻译吧。karyotype = data/karyotype/karyotype.human.txt #核型文件，包括染色体具体信息和颜色配置，

为染色体排序，博主认为这个不是很重要，要记那么多模式，还不累死啊！还不如老老实实地直接写全比较好，当然心有余力怎么着都行，主配置文件有变化：>>>>karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default

过滤这一节主要讲的是有时候不需要全部的染色体都出现，那么就需要过滤一部分了，其他配置见第二节，主配置文件如下：>>>>karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = nochromoso

label的标签有的时候没有办法完全表现出来，这个主要还是由于空间不够造成的，只看想要的基因还是很不错的，主配置文件：karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = nochromosomes =

承接上一篇，只有主配置文件有变化，这次的内容为规则与链接，副内容为贝塞尔曲线：karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = nochromosomes = /hs[1-4]$/chromos

染色体之间的间隔也可以控制，根据上一节文件，来配置主配置文件和spacing文件：>>>>karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = nochromosomes = h

最简单的数据结构的应用，把相同id的行合并，涉及引用等知识：#!/usr/bin/env perluse warnings;use strict;my %hash;open FA, $ARGV[0] || die $!;while (){ chomp; my @tmp = split; my $string = "$tmp[1]\t$tmp[2]\t$tmp[3]\t$tm

WGCNA（Weighted Correlation Network analysis）是一个基于基因表达网络权重构建，描述基因表达的关联模式的R包。挺拗口的吧，其实简单点的话分析基因的共表达网络，就是两个样本有表达量，那么博主根据表达量可以计算相关性，但如果加入一些新的权重，比如重量、高度、应激条件等等，相当于把基因表达与条件结合起来分析两者之间的关联性或相关性，当然表达量是最关键的。由此也可以

这一节将为染色体加上个标签，方便处理染色体的信息，比如移动位置：>>>>karyotype = data/karyotype/karyotype.human.hg19.txtchromosomes_units = 1000000chromosomes_display_default = nochromosomes = hs1[a]:1-50;

好吧，上次由于交代不是很清楚，搞了个抽样放回的程序，今天重新弄了一下，把它写成模块的形式：package Sampling;use warnings;use strict;require Exporter;our @ISA = qw(Exporter);our @EXPORT_OK = qw(swor fac); #swor -- sampling without replacem

博主认为重新写个整体的脚本比较合适，一段一段写容易混淆。层次聚类先将矩阵转置，然后求powers（阈值），根据阈值进行聚类：library(WGCNA)options(stringsAsFactors=FALSE)enableWGCNAThreads()myfile=read.table("da1.nom",sep="\t",header=TRUE)mydata=as.data.fra

#!/usr/bin/env perluse warnings;use strict;die "perl $0 \ne.g. perl $0 96 100 2 > 96.txt \n" unless @ARGV eq 3;my $sample = $ARGV[0]; # 样本数目my $n = $ARGV[1]; # 取出样本的子集数目my $ex_samp = $ARGV

好吧，该到最后的时刻了，识别novel可变剪接的类型，首先要在原理上判断类型的情况，把类型画出来比较好一些：以上即为大致的分类，不过有两点要说明的：intron retained 1用tophat的结果无法实现，这个看覆盖度即可；另外还有两个分类：intergenic和other，顾名思义基因间区域内的和无法识别的，当然也可以按照其他分类，就不多说了，就以此为例吧：#!/usr/bin

想把之前做的可变剪接模型给大家说一下，看看有什么遗漏的没有，由于当时想法比较复杂，所以程序有点多，大致分三个部分来进行。首先，拿到的结果是tophat给出的junction的数据，其次博主使用的数据库是ensembl的数据库，gencode也可以，先得到已知的参考junction：#!/usr/bin/env perluse warnings;use strict;die "per

博主其实对未知的可变剪接分类有些困惑，但想了很久才使用了一项比较复杂的算法，接下来并不是分类，而是先转换数据库，因为ensembl|gencode数据库的格式并不能满足博主的需求，转换后更能方便地处理接下来的工作：#!/usr/bin/env perluse warnings;use strict;my (%gene);open GTF, $ARGV[0] or die $!;wh

可变的半径，至少是辐射状的图像吧，其实比较简单其他配置文件未变，主配置文件如下：>>>>karyotype = data/karyotype/karyotype.human.hg19.txtchromosomes_units = 1000000chromosomes_display_default = yeschromosomes_radius

又是label？？？确实这一节还是label，不过具体讲了一些细节东西，比如字体的改变之类的，一般没人想去改它吧- -！label配置文件和主配置文件的变化：>>>show_ticks* = no #不显示ticks>karyotype = data/karyotype/karyotype.human.txtchromosomes_units

终于到这一章的最后一节了，还是基础讲解。角度偏差在circos的设置以地球方向为例：北-90，南90，东0，西180，当然几乎所有的定位默认是-90。举个例子来说，之前我们写的所有程序几乎都是从chr1开始，有没有发现chr1的开始位置就是-90呢？所以想要从哪个方向开始，就在里面重写吧。其次我们之前的程序都是顺时针clockwise，当然还有逆时针的countclockwise，主配置文件如

呵呵，好久没更新了。。。话说这段时间比较忙，自己用之前的三章画了一些circos图深有感触，算是实战了，自我感觉画的还行，但暂时不能放上来，因为还没发表。。。这一章节是高亮显示数据，显示相应区域的数据或重新定义一个新的ideogram很有用处，结构也比较简单：#bands.confshow_bands = yesfill_bands = y

紧跟上一节的内容，主配置文件有变化，其实就是将输入文件直接格式化成颜色位置等等，没有必要一步一步的画了，文件格式如下：hs1 1298972 1300443 fill_color=bluehs1 1311738 1324571 fill_color=red,r0=0.6r,r1=0.6r+50phs1 1397026 1421444 fill_color=green,r0=1.1r,r1=

z深度，之前已经解释过了，算是优先级，值越大覆盖值小的颜色，当然z深度也有其他的功能，主配置文件如下：>>>>karyotype = data/karyotype/karyotype.human.txtchromosomes_units = 1000000chromosomes_display_default = nochromosomes