您的位置:新葡亰496net > 电脑系统 > 新葡亰496net:基因组干货,的一对脑洞操作

新葡亰496net:基因组干货,的一对脑洞操作

发布时间:2019-10-13 11:55编辑:电脑系统浏览(191)

    把当下文件夹的文本名用","连接成一行,恐怕将多行调换为一行

    ls | paste -s -d ","  # -s 选项将输入进行一次性粘贴
    ls | xargs | sed 's/ /,/g'  #xargs 将输入作为参数(空格分隔)传入
    ls | awk '{printf "%s,",$0}'
    

    本次讲awk 的数组及采用

    申明:本学习笔记指标用于笔者个人学习,其内容整理出自MOOC电影大学数学实验张勇先生等的教程课件,再次申明,请勿转发。www.icourse163.org/learn/UESTC-235004

    明日讲对于基因注释.gtf格式文件举行操作。

    将行逆序输出

    sed '1!G;h;$!d'file  # 1!G 第一行不执行G命令,从第二行开始执行;$!d 最后一行不删除;第一行自动存入模式空间,将模式空间内容(第一行)放到保持空间(h),然后删除模式空间内容(d,否则它会自动输出),第二行自动存入模式空间,(开始用G)将保持空间(第一行内容)接到模式空间(第二行)后,将当前模式空间(第二行 第一行)放到保持空间(h),然后删除当前模式空间(d),依次类推,最后一行不删除模式空间,再自动输出模式空间内容
    tac file
    

    awk 中数组叫做关联数组(associative arrays),因为下标能够是数字也足以是字符串。awk 中的数组不必提前注明,也不必注明大小,直接赋值就行。


    人类基因组注释文件在测序解析中作为那么些主要的组装参照物,並且在其中包蕴了差不离人类基因组全数的着力音信,我们从GENCODE(https://www.gencodegenes.org/) 下载gtf格式的注明文件,

    新葡亰496net 1

    删除#伊始的注释行

    sed '/^#.*/d' test.txt
    
    # 可以用数值作数组索引(下标)
    Tarray[1]=“cheng mo”
    Tarray[2]=“800927”
    
    # 可以用字符串作数组索引(下标)
    Tarray[“first”]=“cheng ”
    Tarray[“last”]=”mo”
    Tarray[“birth”]=”800927”
    # 也可以这样:
    Tarray[“first”t"last"]=”chengmo”
    

    1.1骨干语法

    此地下载hg19版本的gtf文件:

    新葡亰496net 2

    # 使用wget命令下载注释文件:
    wget -c -t 0 ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
    # -c 表示端点续传
    # -t 表示反复尝试的次数,0为不限次数 
    # 解压缩并删除原文件:
    gunzip ./gencode.v19.annotation.gtf.gz
    # 查看这个gtf文件大小,注意最后还有个点,表示当前目录:
    du -a -h -d 1 .
    # -a表示显示当前目录下每个文件的磁盘使用情况
    # -h表示以以K,M,G为单位,提高信息的可读性
    # -d N与--max-depth=N一样,表示只深入到第N层目录,设置为0,即表示不深入到子目录
    

    新葡亰496net 3

    去掉每行带头4个字符

    cut -c 4- test.csv
    
    或许先成立bed文件:
    # 使用 randomBed创建bed文件:
    randomBed -seed 2 -n 1000 -l 400 -g ./hg19.chrom_24.sizes >./a.bed
    # 查看前10行:
    head a.bed
    

    新葡亰496net 4

    1.变量命名准则

    a.必得以字母开始;b.区分轻重缓急写;c.可有字母、下划、线数字组合

    example: numcar or num_car; %变量名应能够反映其实际意义。

    GTF全称为gene transfer format,那一个GTF注释文件收缩后独有40M,但解压出来后有1.1个G,消息及其冗余,内容格式是这么的,GTF文件都由9列数据整合,字段必需tab分开,除了每一行的终极贰个字段之外,都无法不带有有值;空驶列车要求用斯洛伐克共和国(The Slovak Republic)语点号.标示,具体每一列的注解如下:
    1. seqname - 染色体或scaffold名称,染色体名称能够未有"chr"前缀,何况那些名称必需时ensemble中的。
    2. 讲授来源
    3. 队列类型
    4. start - 初始为点,1-base格式。
    5. end - 终止为点,1-base格式.
    6. score - 浮点型数值.未有就用.表示
    7. strand - (forward) 或 - (reverse).
    8. frame - One of '0', '1' or '2'. '0' 表示种类的率先个碱基初叶是密码子, '1'象征系列的第三个碱基起首是密码子,依此类推。未有就用.表示
    9. attribute - 封号分割的连串特征音信,包含基因ID,转录本ID,基因类型等.
      更实际详细的新闻能够查阅这里:https://www.gencodegenes.org/data_format.html

    新葡亰496net 5

    对文本首先列实行总括

    awk -F "," '{count[$1]  } END{for (record in count) print record,count[record] }' test.csv  #count[$1]  创建关联数组count[$1]并进行计数
    
    以此bed文件突显的第四列是行号排序,今后我们想将它替换为相应的染色体的长短,举例chr1的长度是249250621,chr2的长度是243一九九一73,那一个新闻是在hg19.chrom_24.sizes以此文件里有些:
    # 查看hg19.chrom_24.sizes
    cat hg19.chrom_24.sizes
    

    新葡亰496net 6

    2.赋值语句

    变量名=表达式;

    example: a=[2,3,4,5]; or a(2)=10;

    其他赋值语句形如:变量名=函数名(输入参数列表);[变量名列表]=函数名(输入参数列表)

    在这里个gtf文件中隐含人类基因组全体的protein-coding, lncRubiconNA, mi智跑NA, pseudogene,sno揽胜极光NA等富有转录本消息,gtf文件的9列之间都是Tab分割的,但结尾一列特征列,内容之间是用封号空格细分的.awk 暗中同意分割是[t ] 也正是三个或多个Tab或空格,这里大家为了更加好的对GTF文件的第九列分割,重新定义FS(田野 separator,分割符)为"t|(; )",表示Tab或封号空格 分割。NENVISION表示Number of row,也正是行号,如下图,GTF文件前5行是无需的,所以加个判断规范NR>5

    新葡亰496net 7

    上边包车型大巴通令提取了第1列的染色体, 第9列的geng_id和第10列的转录本ID音信:

    awk 'BEGIN{
        FS="t|(; )"
        }
        {
        if(NR>5){
            print $1,$9,$10
            }
        }' gencode.v19.annotation.gtf | head 
    

    新葡亰496net 8

    遵照地点定义的分割符,1-8列的新闻还和下面的一律,但9-13的新闻已然是gene_id , transcript_id , gene_type , gene_status , gene_name就这样类推。

    如此这般我们就能够提取想要的信息了,比方笔者只想要蛋白编码基因的新闻,但不包蕴转录本:

    awk 'BEGIN{
        FS="t|(; )";
        #设置输出分割符为Tab,默认为空格
        OFS="t"
        }
        {
        #注意要加反斜杠对双引号转义:
        if(NR>5 && $3=="gene" && $11=="gene_type "protein_coding""){
            print $1,$4,$5,$13,$7,$11
            }
        }' gencode.v19.annotation.gtf >./protein_coding_gene.txt
    # 使用gedit 查看一下:
    gedit ./protein_coding_gene.txt
    

    新葡亰496net 9

    诸有此类大家就获得了讲明文件中的全数203五十多个蛋白编码基因的地方消息及基因名。但大家还要求越来越对于出口音信举办优化,举例输出时第四列去掉"gene_name"那样的字段,只保留后边的称谓,大家下一次勇往直前。
    更加多原创精粹内容敬请关心生信故事集

    新葡亰496net 10

    对文件第四列用":"切割成两列并将最后一列结果 1,然后输出全体列

    awk -F "," '{split($4,array,":");print $1,$2,$3,array[1],array[2] 1}' test.csv  #split切割$4存到数组array中,array[1]和arrya[2]即为切割后的两个区域
    
    轮换后的结果应当是以此样子的:

    新葡亰496net 11

    3.表明式语句

    二个讲话能够唯有表明式,系统自动将表明式的结果赋值给MATLAB内部变量"ans"。

    对文本第二列求均值

    awk -F "," '{sum =$2} END {print "Average = ", sum/NR}' test.csv
    
    第四列已经替换为了呼应染色体长度。

    4.语句分隔符

    支行和逗号,要是不加分号,系统交易会示运算后的结果。

    兑现DNA体系反向互补

    cat seq.txt | sed 'y/ATGC/TACG/' |rev
    
    代码及注释如下:
    awk 'BEGIN{OFS="t"}{
        # 在按行读第一个文件时将第二列数据存储到以第一列染色体名为下标的数组a里:
        if(NR==FNR){
            a[$1]=$2
            }
        # 读第二个文件,注意其中的a[$1]已经是a.bed中的$1.
        if(NR>FNR){
            print $1,$2,$3,a[$1],$5,$6
            }
        }' hg19.chrom_24.sizes a.bed >b.bed
    

    5.常用命令、快速键

    clear          清除工作空间中的变量

    如:clear   变量名列表

    示例:

    clear A B    清除变量A,B

    clc              清除命令窗口内容

    who            列出当前工作空间有所变量名称

    whos          列出当前工作空间变量更加的多音信(维数,占用内部存款和储蓄器字节数等)

    whos          变量名列表

    示例:

    whos v1 v2列出变量v1,v2的越来越多音信

    连忙键:向上方向键、向下方向键

    用以浏览命令窗口历史命令、语句


    某一行插入别的一个文本的开始和结果

    sed '2 r a.txt' test.csv 
    
    仍是能够透过length函数获得数老总度,但awk有几许个本子,用的可比多的时gawk,在使用函数以致awk内置变量时应有安转gawk:
    sudo apt-get install gawk
    awk '{a[$1];print "当前数组长度为:",length(a)}' hg19.chrom_24.sizes
    

    新葡亰496net 12

    乘势读取行数大增,数首席营业官度也在大增。


    越多原创雅观内容敬请关怀生信散文

    新葡亰496net 13

    1.2数组的创制与使用

    对二个文件遵照第一列进行筛选,筛选规范是必需在别的二个文本的第一列出现过

    awk -F "," '{if(NR==FNR){count[$1]=1}else if(count[$1]==1){print $0}}' chr.txt test.csv  #将第一个文件第一列的值存入关联数组,并给值为1,如果第二个文件建立的关联数组对应值为1,说明在第一个文件第一列出现过,则输出整行
    

    1.2.1创办数组

    1. 选用方括号

    同一行的因素用“空格或逗号”分隔,不一致行的要素用“分号或换行”分隔。

    如: X=[1,2,3;4,5,6;7,8,9];

            Y=[2,3,4

                  5,6,7];

    2.冒号符操作

    用来创制行向量a:step:b 当中a:b等同于a:1:b,a为初始值,step为增量,b用于判断向量终点值。

    x=1:5 表示x=[1 2 3 4 5],增量默感觉1

    x=1:2:9 表示x=[1 3 5 7 9]

    x=10:-2:1 表示x=[10 8 6 4 2]

    3.linspace(a,b,n)

    n-1等分区间[a, b]的节点组成的行向量(含区间端点a, b)

    示例:x=linspace(-2, 2, 5) %表示x=[-2 -1 0 1 2]

    万一要爆发叁个间距上的均匀节点,况且内定所发出数组的成分个数,则应用linspace更为有扶持。

    4.拼接

    演示格式1:[A B] 横向拼接供给A,B行数一样,

    亲自过问格式2:[A; B] 纵向拼接,需求A,B列数同样.

    示例:z=[rand(2,3), rand(2,2) ]

    5.空矩阵[ ] 产生贰个空矩阵

    示例:a=[ ]

    6.调用函数创设

    a = zeros(m, n) 产生叁个m行、n列的零矩阵;多用于变量的最早化

    a = ones(m, n) 发生二个m行、n列的因素全为1的矩阵

    a = eye(m, n) 发生四个m行、n列的单位矩阵

    对文本第二列和第三列实行举办

    开展前四列
    新葡亰496net 14
    扩充后改为三列
    新葡亰496net 15

    awk -F "," '{for (i=$2;i<=$3;i  ) {print $1,i,$4}}'  test.csv 
    

    1.2.2修改和领取数组中的成分

    由此数组下标访谈:(1)下标为大于等于1 的整数;(2)下标不能够越界

    常用语法:示例:x(i), x(a : b : c), x([a b c d]),x(i,j)

    获取子阵:

    猎取某一行 A(r, :) 第r行;

    获得某一列 A(:, c) 第c列

    收获子阵A(行下标集,列下标集)

    修改成分:用赋值语句修改。要是赋值语句右边不是二个标量,则须要赋值语句两边表示的数组维数要长久以来,不然变成维数不等同的不当。

    参考用法1:A(i,:)=B(k,:),

    参谋用法2:A([1 2],:)=V

    对八个公文相继merge

      这里四个文件行数相等,个中ampl列将新的和旧的染色体、地方关系起来,第一个公文将第五列(ampl列,值为ampl1,ampl2...)存入一二三列(旧染色体,旧起头地方,旧甘休地点)为下标的关联数组ampl,第1个文本根据一二三列(旧染色体,旧起头地方,旧停止地点)抽出关联数组的值(ampl1,ampl2...),将涉及数组的值作为关周详组下标新创制关联数组Ampl,将第叁个文本的值(1,2,3,4,5列,当中4、5列是大家要的音讯)用sprintf生成字符串存入Ampl,第三文书依照第四列(ampl1,ampl2...),用split切割sprintf生成的字符串,抽取第三个公文存入的值(这里只收取了须求的4,5列,123列的值输出第多个文本的123列(新染色体,新开第4地点,新竣事地方)的值)。那样Oldpanel_start_end.sort.bed 对应的旧的染色体和职分,被hg38amplicon_start_end.bed新的叁个染色体和岗位替代,而且将旧文件染色体和任务在amplGChg19.txt 对应的新闻成功转移到新生成的新岗位文件中

    awk 'BEGIN{FS="t";OFS="t"}{if(NR==FNR){ampl[$1,$2,$3]=$5;N=NR}else if(NR<=2*N){Ampl[ampl[$1,$2,$3]]=sprintf("%s,%d,%d,%s,%s",$1,$2,$3,$4,$5);}else{split(Ampl[$4],array,",");print $1,$2,$3,array[4],array[5],$4}}' Oldpanel_start_end.sort.bed amplGChg19.txt hg38amplicon_start_end.bed | sort -k1 > hg38amplicon_Gene_GC.txt
    

    1.2.3删减数组中的成分

    操作办法:将空矩阵赋值给相应子阵达到删除指标。

    用法:

    A(i1:i2,:)=[]%删除A由i1:i2钦命的行

    A(:,j1:j2)=[]%删除A由j1:j2钦命的列

    对多少个文件去重取并集

    cat NewpanelGene.bed Oldpanel.gene.bed | sort -u > merge.gene.bed  #sort -u = sort | uniq ,相当于sort 之后,将重复相邻行变成只有一行
    

    1.2.4 end在存取数组成分方面包车型大巴非正规用法

    用法:end在下标表明式中代表最后三个下标值

    假诺end出现在三个向量的下标中,则表示向量的要素个数。

    假诺end出现在二个矩阵的行下标地点,则意味矩阵的行数。

    倘若end出现在多个矩阵的列下标位置,则意味着矩阵的列数。

    示例:x=[1 5 9; 2 6 10; 3 7 11; 4 8 12];

    x(end,2)= 0; x%将矩阵x的末梢一行第2列成分赋值为0

    任何例子:

    -------------------操作向量示例

    t = rand(1,10);

    x1 = t(1:end-1) %取第1个-倒数第2个

    x2 = t(end-2:end) %取尾数第三个-尾数第4个

    –操作矩阵示例

    A = rand(3)

    B = A(1:end-1, : ) %取A的第1行-倒数第2行

    C = A(:, [2:end]) %取A的第2列-倒数第1列


    对文件依照标识起先的行开展分割

    比如
    新葡亰496net 16

     awk '/>chr/{split($0,array,">");out=array[2]};{print > out}' test.fa
    

    出口chr1,chr2多个文本

    1.3运算符

    出口文件奇数行和偶数行

    sed -n 'p;n' test.txt #输出奇数行
    sed -n 'n;p' test.txt #输出偶数行
    

    1.3.1算术运算符

    矩阵转置B. '矩阵共轭转置B'

    矩阵加减:A B,A-B,A与B维数一样或内部之一为标量

    矩阵相乘:A*B,A与B为矩阵或内部之一为标量

    矩阵左除:AB,当A为方阵表示: A-1B

    矩阵右除:A/B,当B为方阵表示AB-1,或B为标量

    矩阵幂:A^n,A为方阵

    数组对应成分总结:

    数组相乘:C=A.*B

    数组右除:C=A./B;

    数组左除:C=A.B

    数组幂:C=A.^B

    渴求:A, B同维数或内部之一为标量

    新葡亰496net:基因组干货,的一对脑洞操作。统计GC含量

    echo "TTCCTTGAAATAAGTGTGATT" | awk '{s=gsub("[GC]","N",$0);print s/length}'
    

    1.3.2涉及运算符

    新葡亰496net 17

    关联运算符表明

    去除windows换行符

    cat test.txt | sed 's/r//g' 
    

    1.3.3逻辑运算符

    逻辑运算的值为0(代表“假”)或1(代表“真”)

    三种运算符:

    与(and) &

    或(or) |

    非(not) ~


    1.4变量数据类型

    器重的数据类型:double char sym struct cell

    a=rand(3); b='Li San';%a为double型,b为char型

    syms x, y=1 x^2 %x,y为sym类型;对y赋值的语句含符号对象

    F.name='li San', F.birth=1999, F.src=rand(3)%F为struct型

    whos a b x y F

    Name Size Bytes Class Attributes

    F 1x1 620 struct

    a 3x3 72 double

    b 1x6 12 char

    x 1x1 112 sym

    y 1x1 112 sym

    查阅变量类型

    示例:

    a=rand(3); b='abc'

    class(a), class(b)

    新葡亰496net:基因组干货,的一对脑洞操作。运营结果:

    ans= double

    ans= char

    运用函数class

    用法:class(变量名)

    该函数重返变量的数组类型的char型数组,如'double', 'char'。

    cell数组基本用法

    创建数组用法:

    a=cell(m,n)

    存取cell数组用法示例:

    a{i} i为下标

    a{i,j} i,j分别为行、列下标

    特点:壹个cell数组中的元素的花色能够互差异


    1.5主导输入与格式化输出操作函数

    八个函数:

    input 输入函数

    输入函数input

    先是种用法:

    input(提醒消息字符数组)

    用以输入通常品种数据

    第三种用法:

    input(提示字符串,'s')

    用以输入字符数组(含首个参数's')

    disp呈现数组内容函数

    来得数组的要素

    数组呈现函数disp(变量名)

    特点:展现数组内容,但不出口变量名

    多用来调节和测量检验程序时呈现数组内容

    sprintf将数组内容格式化为字符串

    格式化输出函数sprintf

    效率:将数据格式化输出为字符串

    用法:str = sprintf(formatSpec,A1,A2,...,An)

    将数组A1,A2,...,An根据参数formatSpec格式化为字符串赋给str.

    %d 格式化整数%f 格式化浮点数

    %c 格式化单个字符%s 格式化字符数组

    百分号符号字符后能够加个整数, 用以限制输出化为字符串的尺寸,比方:], %5s。


    本文由新葡亰496net发布于电脑系统,转载请注明出处:新葡亰496net:基因组干货,的一对脑洞操作

    关键词:

上一篇:新葡亰496net:一声令下大全

下一篇:没有了