精华内容
下载资源
问答
  • Usage: seqkit subseq [flags] Flags: --bed string by tab-...#根据bed、gtf文件提取基因 seqkit subseq --bed bedfile.bed -o gene.fa genomefile.fa seqkit subseq --gtf gtffile.bed -o gene.fa genomefile.fa
  • fasta序列操作神器——seqkit

    千次阅读 2019-10-04 20:48:57
    seqkit seq test.fa -r > test_re.fa 2.取互补序列 seq test.fa -p > test_com.fa 3.取反向互补序列 seqkit seq test.fa -r -p > test_re_com.fa 4.DNA序列转换为RNA序列 seqkit seq test.fa --nda2rna &.....

    一、序列操作:

    1.取反向序列

    seqkit seq test.fa -r > test_re.fa

    2.取互补序列

    seq test.fa -p > test_com.fa

    3.取反向互补序列

    seqkit seq test.fa -r -p > test_re_com.fa

    4.DNA序列转换为RNA序列

    seqkit seq test.fa --nda2rna > test_rna.fa

    5.RNA序列转换为DNA序列

    seqkit seq test.fa rna2dna > test_dna.fa

    6.将序列以小写字母的形式输出

    seqkit seq test.fa -l > test_lower.fa

    7.将序列以大写字母的形式输出

    seqkit seq test.fa -u > test_upper.fa

    8.指定每行序列的输出长度(为0的话,代表为一整行,默认的输出 长度是60个碱基)

    seqkit seq test.fa -w 10 > test_10.fa (指定序列的长度为10)

    9.将多行序列转换为一行序列

    seqkit seq test.fa -w 0 > test_w.fa

    10.只输出序列

    seqkit seq test.fa -s -w 0 > test_seq.fa

    11.将只输出的序列的,指定每行输出的碱基数

    seqkit seq test_seq.fa -s -w 40 > test_seq40.fa

    注意10,11的微妙之处

    11,12也可以一步完成:

    seqkit seq test.fa -s -w 20 -o test_20.fa

    二、Fasta/q之间以及与tab格式互换

    10.将fataq文件转化为fasta格式.

    seqkit seq fq2fa test.fq -o test.fa

    11.将fasta格式转化为tab格式

    seqkit fx2tab test.fa > test_tab.fa (没有seq参数)

    三、序列信息统计

    1.序列碱基含量

    seqkit fx2tab -l -g -n -i -H test..fa (这些参数组合起来比较好看)

    2.序列长度的整体分布统计

    seqkit stat test.fa
    seqkit grep [flags]

    参数:

    -n, --by-name

    匹配整个序列的名字,包含deion部分,而不是序列id。

    -s, --by-seq

    匹配序列

    -d, --degenerate

    pattern/motif 包含简并碱基

    -i, --ignore-case

    忽略大小写

    -v, --invert-match

    输出不匹配此模式的内容

    -p,

    匹配模式,支持连续写多个模式,匹配任一模式即输出。如-p ^ATG -p TAA$。注意该功能仅能正向匹配,不能实现对互补链匹配。

    -f, --pattern-file string

    支持匹配模式写到一个文件中,如要提取的序列ID。

    -R, --region string

    匹配位置选择。e.g 1:12 for first 12 bases, -12:-1 for last 12 bases

    -r, --use-regexp

    使用正则表达式,必须加入此参数,如^匹配首端。同-p联合使用。

    举例:

    seqkit grep -s -r -i -p ^atg cds.fa#选取有起始密码子的序列

    seqkit grep -f list test.fa > new.fa#根据ID提取序列

    seqkit grep -s -d -i -p TTSAA#简并碱基使用。S 代表C or G.

    seqkit grep -s -R 1:30 -i -r -p GCTGG##匹配限定到某区域

    五、motif定位

    对grep的拓展,可以正反链同时匹配,输出匹配的位置。

    seqkit locate [flags]

    参数:

    -d, --degenerate

    pattern/motif contains degenerate base

    -i, --ignore-case

    ignore case

    -P, --only-positive-strand

    only search at positive strand

    -p, --pattern value

    search pattern/motif

    -f, --pattern-file string

    pattern/motif file (FASTA format)

    举例:

    seqkit locate -i -d -p AUGGACUN test.fa

    输出结果:

    seqID

    patternName

    pattern

    strand

    start

    end

    matched

    cel-mir-58a

    AUGGACUN

    AUGGACUN

    81

    88

    AUGGACUG

    ath-MIR163

    AUGGACUN

    AUGGACUN

    122

    129

    AUGGACUC

    六、多个序列文件比较寻找相同的序列或者ID相同的序列

    seqkit common [flags]

    参数:

    -n, --by-name

    匹配整个序列的名字,包含deion部分,而不是序列id

    -s, --by-seq

    match by sequence

    -i, --ignore-case

    ignore case

    -m, --md5

    use MD5 reduce memory usage

    举例:

    1、By ID (default,>后面,空格之前的名字)输出ID名字相同的。

    seqkit common test1.fa test2.fa -o common.fasta

    2、By full name(整个序列的名字,包含deion部分)。输出序列名字相同的。

    seqkit common test1.fa test2.fa -n -o common.fasta

    3、输出要比较的文件中序列相同的序列

    seqkit common test1.fa test2.fa -s -i -o common.fasta

    4、输出要比较的文件中序列相同的序列 (for large sequences)

    seqkit common test1.fa test2.fa -s -i -o common.fasta --md5

    七、提取部分序列

    如随机抽取10000条FASTQ序列做NT污染评估。同时他也可以对FASTA序列提取

    seqkit sample [flags]

    参数:

    -n, --number int

    sample by number (result may not exactly match)

    -p, --proportion float

    sample by proportion(按比例提)

    -s, --rand-seed int

    rand seed for shuffle (default 11)

    -2, --two-pass

    2-pass modelower memory

    举例:随机抽取序列

    seqkit sample -n 10000 -s 11 test1_1.fq -o sample.fq

    seqkit sample -p 0.1 -s 11 test1_1.fq -o sample.fq

    八、排序输出命令

    seqkit sort [flags]

    参数:

    -l, --by-length

    按照序列长度排序

    -n, --by-name

    by full name

    -s, --by-seq

    按照序列排序

    -i, --ignore-case

    按序列排序时忽略大小写

    -r, --reverse

    反向排序

    -2, --two-pass

    对于FASTA序列排序可以减少内存

    举例:

    seqkit sort -ltest.fa

    九、文件切割

    seqkit split [flags]

    参数:

    -i, --by-id

    split squences according to sequence ID

    -p, --by-part int

    将一个文件分割成N 份

    -s, --by-size int

    将一个文件按照N 条序列一个文件进行分割

    -O, --out-dir string

    output directory (default value is infile.split)

    -2, --two-pass

    two-pass mode to lower memory usage(only FAST)

    举例:

    seqkit split hairpin.fa.gz -p 4

    转载于:https://www.cnblogs.com/huangyinger/p/10421805.html

    展开全文
  • 1.orthofinder介绍OrthoFinder是一种快速、准确和全面的比较基因组学分析工具。它可以找到直系和正群,为所有的正群推断基因树,并为所分析的物种推断一个有根的物种树。OrthoFinder还为比较基因组分析提供全面的...

    1.orthofinder介绍

    OrthoFinder是一种快速、准确和全面的比较基因组学分析工具。它可以找到直系和正群,为所有的正群推断基因树,并为所分析的物种推断一个有根的物种树。OrthoFinder还为比较基因组分析提供全面的统计数据。OrthoFinder使用简单,只需运行一组FASTA格式的蛋白质序列文件(每个物种一个)。

    2.基础知识介绍

    Orthologue(直系同源基因)指的是来自两个物种的基因。Orthologue是由两个物种的最后共同祖先(LCA)上的单个基因进化而来的成对基因(图1A和B)。正群是同源概念在物种群中的自然延伸。一个Orthogroup(正交群)是由一个物种的LCA中的单个基因进化而来的一组基因(图1A)。当观察基因树时,一个邻位群体中基因的第一次分化是一个物种形成事件,对同源基因来说也是如此。

    作为基因复制事件的结果,当观察直系同源基因和正交群时,可能会有来自同一物种的多个基因。在这个例子中(图1A和B),人类和老鼠的HuA基因是鸡中ChA1和ChA2的同源基因。再看一下正交群,我们发现有两个鸡的基因(图1A),但是只有一个来自老鼠和人类的基因。一些作者将ChA1和ChA2基因作为HuA的共同源基因,以强调存在多个同源基因的事实。由于基因重复和丢失在进化中经常发生,一对一的直系同源物很少见,通过分析正交群所有直系同源的情况(一对一,多对一,多对多),我们可以分析数据的所有情况。

    paralogues (旁系同源基因)是指在基因复制事件中从单个基因中分离出来的成对基因,鸡的两个基因ChA1和ChA2是旁系同源基因(图1C)。来自不同物种的两个基因如果在基因重复事件中彼此分离,也可能是同源的。由于基因树中所有的分支事件要么是物种形成事件(产生直系同源基因),要么是重复事件(产生旁系同源基因),因此同一正交群中任何不是直系同源基因的基因必然是旁系同源基因。

    52c2b99615f6

    图1

    直系同源物是同源性基因,是物种形成事件的结果。Paralogs(旁系同源物)是同源基因,是重复事件的结果。可以看到(图2),不同物种间的α-chain gene互为Orthologs(直系同源物)。正交群用来形容自一组物种的LCA中的单个基因的基因组(α-chain gene)。然后同一物种间α 和β chain gene互为Paralogs(旁系同源物)。最后所有这些关系都可以由OrthoFinder来识别。

    52c2b99615f6

    图2

    3.安装

    在这里推荐大家使用conda来安装,简单明了,不用担心其他Dependencies的安装

    $conda install orthofinder

    4.运行orthofinder

    运行Orthofinder,相当简单的操作, "-f"输入目录,里面包含你需要运行的蛋白质fasta文件,“-t”所用到的CPU数目,“-S”选择的比对模式(默认diamod,可选blast)。

    我使用的物种有九个:大豆(G.max)、菠萝(A.comosus)、柑橘(C.sinensis)、无油樟(A.trichopoda)、拟南芥(A.thaliana)、黄瓜(C.sativus)、水稻(O.sativa)、小立碗藓(P.patens)、江南卷柏(S.moellendorffii)

    将所有物种的蛋白文件放到Angiospermae文件夹下

    $nohup orthofinder -t 16 -f Angiospermae/ -S blast &

    生成的结果会存储于Orthofinder/Results_XXX文件中,现在简单看看里面有啥。

    Citation.txt Orthologues

    Comparative_Genomics_Statistics Phylogenetically_Misplaced_Genes

    Gene_Duplication_Events Putative_Xenologs

    Gene_Trees Resolved_Gene_Trees

    Log.txt Single_Copy_Orthologue_Sequences

    Orthogroups Species_Tree

    Orthogroup_Sequences WorkingDirectory

    我们主要使用Orthoogroups查看正交群的基因和使用 Single_Copy_Orthologue_Sequences里的单拷贝基因构建系统发育树

    WorkingDirectory其中包含运算过程的中间文件,例如blast结果,如果我们想去掉某一物种,在SpeciesIDs.txt中将该物种注释掉

    $vi SpeciesIDs.txt

    #0: Acomosus.pro.fa

    1: Athaliana.pro.fa

    2: Atrichopoda.pro.fa

    3: Csativus.pro.fa

    4: Csinensis.pro.fa

    5: Gmax.pro.fa

    6: Osativa.pro.fa

    7: Ppatens.pro.fa

    8: Smoellendorffii.pro.fa

    ~

    $nohup orthofinder -b WorkingDirectory

    如果我们想增加额外的物种进行分析,可以按照如下方式运行

    $nohup orthofinder -b WorkingDirectory -f new_fasta_directory

    5.利用单拷贝基因构建系统发育树

    Orthofinder的 Single_Copy_Orthologue_Sequences下存放着单拷贝同源基因的序列,我们可以利用这些序列构建系统发育树

    5.1.多序列比对

    多序列比对推荐使用muscle

    $wget http://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz

    $tar xzvf muscle3.8.31_i86linux64.tar.gz

    $chmod +x muscle

    将muscle添加至环境变量

    单拷贝的序列较多,所以使用shell批量运行

    $vi bash.sh

    for i in *.fa

    do muscle -in $i -out $i.1

    done

    $sh bash.sh

    5.2.提取保守序列

    *.1文件是比对好的序列文件,接下来使用Gblocks提取保守序列

    $wget http://molevol.cmima.csic.es/castresana/Gblocks/Gblocks_Linux64_0.91b.tar.Z

    $tar Zxf Gblocks_Linux64_0.91b.tar.Z

    添加至环境变量

    $vi bash2.sh

    for i in *.1

    do Gblocks $i -b4=5 -b5=h -t=p -e=.2

    done

    $sh bash2.sh

    -t= default:p设置序列的类型,可选的值是 p,d,c 分别代表 protein, DNA, Codons 。

    -b1= default:( 序列条数的 50% + 1 ),设定保守性位点必须有 >= 该值的序列数。该参数后接一个 integer 数,默认下比序列条数的 50% 大 1.

    -b2= default: 序列条数的 85%,确定保守位点的侧翼位点时,其位点必须有 >= 该值的序列数。该值必须要比 -b1 的值要大。

    -b3= default: 8,最大连续非保守位点的长度。

    -b4= default: 10,保守位点区块的最小长度。该值必须 >=2 。

    -b5= default: n, 设置允许含有 Gap 位点。可选的值有 n,h,a 分别代表 None, With Half, All 。 当为 h 时,表示

    -e= default: .2, 设置输出结果的后缀。

    5.3.序列合并

    比对好之后我们需要将所有文件合并,在合并之前,需要对每个序列文件进行排序,并将多行序列转换为一行序列,此时用到的工具是seqkit

    $conda install seqkit

    $vi bash3.sh

    for i in *.2

    do seqkit sort $i >$i.3

    seqkit seq $i.3 -w 0 > $i.3.4

    done

    $sh bash3.sh

    $mkdir new

    $mv *.4 new/

    $cd new

    $paste -d " " *.4 > all.fa

    $sed -i "s\ \\g" all.fa

    此时我们已经把所有的单拷贝序列合并,接下来使用notepad++把所有的ID更改

    52c2b99615f6

    修改前

    52c2b99615f6

    修改后

    5.4.选择合适的替代模型

    修改结束后选择合适的氨基酸替代模型准备构建系统发育树,在这里我使用prottest预测合适的模型,prottest使用phy格式文件,所以用python脚本先将fa文件转换为phy文件,

    import re

    with open('all.fa', 'r') as fin:

    sequences = [(m.group(1), ''.join(m.group(2).split()))

    for m in re.finditer(r'(?m)^>([^ \n]+)[^\n]*([^>]*)', fin.read())]

    with open('all.phy', 'w') as fout:

    fout.write('%d %d\n' % (len(sequences), len(sequences[0][1])))

    for item in sequences:

    fout.write('%-20s %s\n' % item)

    从官网下载prottest

    $ tar zxf prottest-3.4-20140123.tar.gz

    $java -jar /opt/biosoft/prottest-3.4-20140123/prottest-3.4.jar -i all.phy -all-distributions -F -AIC -BIC -tc 0.5 -threads 24 -o prottest.out

    在prottest.out中可以看到最佳模型

    52c2b99615f6

    prottest.out

    5.5.构建系统发育树

    使用raxml构建系统发育树

    $wget https://github.com/stamatak/standard-RAxML/archive/master.zip

    $unzip master.zip

    $cd standard-RAxML

    $make -f Makefile.SSE3.gcc

    $make -f Makefile.SSE3.PTHREADS.gcc

    $make -f Makefile.SSE3.MPI.gcc

    $make -f Makefile.SSE3.HYBRID.gcc

    添加至环境变量

    选择JTT+I+G+F模型构建发育树,外群选择江南卷柏和小立碗藓

    $nohup raxmlHPC-PTHREADS-SSE3 -T 16 -f a -x 123 -p 123 -N 1000 -m PROTGAMMAIJTTF -k -O -o Smoellendorffii,Ppatens /

    -n all.tre -s all.fa &

    运行结束后使用figtree打开RAxML_bipartitions.all.tre,对进化树进行修改。

    52c2b99615f6

    tree

    至此,利用单拷贝基因构建系统进化已经完成

    参考链接

    展开全文
  • 安装软件: 参考「基因组学」使用OrthoFinder进行直系同源基因分析 安装OrthoFinder,然后再安装CAFE wget https://github.com/hahnlab/CAFE/releases/download/v4.2/CAFE-4.2.tar.gz cd CAFE ./configure make ...

    环境准备

    安装软件: 参考「基因组学」使用OrthoFinder进行直系同源基因分析 安装OrthoFinder,然后再安装CAFE

    wget https://github.com/hahnlab/CAFE/releases/download/v4.2/CAFE-4.2.tar.gz
    cd CAFE
    ./configure 
    make 
    make install
    # 如果没有root权限
    mkdir -p ~/opt/biosoft/CAFE-4.2/bin
    make install prefix=~/opt/biosoft/CAFE-4.2
    

    数据准备: 一共会分析mouse, rat, cow, horse, cat, marmoset,macaque, gibbon, baboon, orangutan, chimpanzee, 和 human 12个物种的蛋白数据。既可以在Ensemble上下载,也可以从 https://share.weiyun.com/5ZIjBg8 (密码:3jzdpm)下载。

    解压缩

    tar xf twelve_spp_proteins.tar.gz
    for i in `ls -1 twelve_spp_proteins/*.tar.gz`; do tar xf $i -C twelve_spp_proteins; done
    rm twelve_spp_proteins/*.tar.gz
    

    识别基因家族

    识别物种内和物种间的基因家族分为如下四步:

    1. 仅保留每个基因中有代表性的转录本,去除可变剪切和冗余基因
    2. 建立BLAST数据库,使用blastp进行 all-by-all 的比对
    3. 使用MCL基于blastp结果进行聚类,基因序列相似的通常是一个基因家族
    4. 解析MCL的输出结果,用作CAFE的输入

    将所有最长的转录本合并成单个文件

    提取每个基因中最长的转录本,然后合并成单个文件。

    python python_scripts/cafetutorial_longest_iso.py -d twelve_spp_proteins/
    cat twelve_spp_proteins/longest_*.fa | seqkit rmdup - > makeblastdb_input.fa
    

    All-by-all BLAST

    先用makeblastdb建立BLAST数据库

    makeblastdb -in makeblastdb_input.fa -dbtype prot -out blastdb
    

    之后用blastp进行序列搜索,得到每个序列的相似序列

    blastp -num_threads 20 -db blastdb -query makeblastdb_input.fa -outfmt 7 -seg yes > blast_output.txt &
    

    -seg参数过滤低复杂度的序列(即氨基酸编码为X), -num_threads线程数,此处设置为20.

    这一步运行的时间和你服务器性能有关,你可以读完后续的内容,再继续分析。

    使用MCL进行序列聚类

    现在,我们需要根据BLAST输出中序列相似性信息寻找聚类。这些聚类将是后续用于CAFE分析的基因家族。聚类这一步将通过mcl处理。

    使用shell命令将BLAST转成MCL能够识别的ABC格式

    grep -v "#"  blast_output.txt | cut -f 1,2,11 > blast_output.abc
    

    下一步,创建网络文件(.mci)和字典文件(.tab),然后进行聚类

    mcxload -abc blast_output.abc --stream-mirror --stream-neg-log10 -stream-tf 'ceil(200)' -o blast_output.mci --write-tab blast_output.tab &
    # --stream-mirror: 为输入创建镜像,即每一个X-Y都有一个Y-X
    # --stream-neg-log10: 对输入的数值做-log10 转换
    # -stream-tf: 对输入的数值进行一元函数转换,这里用到的是ceil(200)
    

    根据mci文件进行聚类, 其中主要调整的参数是-I, 它决定了聚类的粒度,值越小那么聚类密度越大,这个值没有想象中的那么至关重要。一般设置为3,你也可以尝试用其他值,然后比较结果。最终的目的是正确分析物种间的直系同源基因。

    mcl blast_output.mci -I 3 
    mcxdump -icl out.blast_output.mci.I30 -tabr blast_output.tab -o dump.blast_output.mci.I30
    

    整理MCL的输出结果

    上一步MCL的输出还不能直接用于CAFE,还需要对其进行解析并过滤。

    第一步是将原来的mcl格式转成CAFE能用的格式

    python python_scripts/cafetutorial_mcl2rawcafe.py -i dump.blast_output.mci.I30 -o unfiltered_cafe_input.txt -sp "ENSG00 ENSPTR ENSPPY ENSPAN ENSNLE ENSMMU ENSCJA ENSRNO ENSMUS ENSFCA ENSECA ENSBTA"
    

    这里的"ENSG00" 是ENSEMBL编号中物种的标识符。并且标识符之前只能有一个空格,否则输出结果就是下面这个情况

    2013053-35aa0a245c33a1d6.png
    错误示范: 输出

    正确的输出结果应该是如下所示

    2013053-ed28373bff474366.png
    正确示范

    第二步,将那些基因拷贝数变异特别大的基因家族剔除掉,因为它会造成参数预测出错。下面的脚本是过滤掉一个或多个物种有超过100个基因拷贝的基因家族,虽然不是特别的严格,但效果和根据拷贝数变异过滤类似

    python python_scripts/cafetutorial_clade_and_size_filter.py -i unfiltered_cafe_input.txt -o filtered_cafe_input.txt -s
    

    将原本的编号替换成有意义的物种名

    sed  -i -e 's/ENSPAN/baboon/' -e 's/ENSFCA/cat/' -e 's/ENSBTA/cow/' -e 's/ENSNLE/gibbon/' -e 's/ENSECA/horse/' -e 's/ENSG00/human/' -e 's/ENSMMU/macaque/' -e 's/ENSCJA/marmoset/' -e 's/ENSMUS/mouse/' -e 's/ENSPPY/orang/' -e 's/ENSRNO/rat/' -e 's/ENSPTR/chimp/' filtered_cafe_input.txt
    sed  -i -e 's/ENSPAN/baboon/' -e 's/ENSFCA/cat/' -e 's/ENSBTA/cow/' -e 's/ENSNLE/gibbon/' -e 's/ENSECA/horse/' -e 's/ENSG00/human/' -e 's/ENSMMU/macaque/' -e 's/ENSCJA/marmoset/' -e 's/ENSMUS/mouse/' -e 's/ENSPPY/orang/' -e 's/ENSRNO/rat/' -e 's/ENSPTR/chimp/' large_filtered_cafe_input.txt
    

    物种树推断

    构建物种树主要分为多序列联配和系统发育树推测两步, 之后在已有进化树的基础上构建超度量树用作CAFE输入。

    多序列联配一般用的是单拷贝的直系同源基因,分别进行多序列联配之后然后合并成单个文件。接着用系统发育树推测软件进行建树,可选软件有

    • 极大似然法: RAxML, PhyML, FastTree
    • 贝叶斯法: MrBayes

    这里不展示具体过程,直接用已有的极大似然树的结果(NEWICK格式),保存为twelve_spp_raxml_cat_tree_midpoint_rooted.txt

    (((cow :0.09289 ,( cat :0.07151 , horse :0.05727) :0.00398) :0.02355 ,(((( orang:0.01034 ,( chimp :0.00440 , human :0.00396) :0.00587) :0.00184 , gibbon:0.01331) :0.00573 ,( macaque :0.00443 , baboon :0.00422) :0.01431):0.01097 , marmoset :0.03886) :0.04239) :0.03383 ,( rat :0.04110 , mouse:0.03854) :0.10918);
    
    2013053-1b4a2857322f754a.png
    已有进化树

    推断超度量树

    超度量树(ultrametric tree)也叫时间树,就是将系统发育树的标度改成时间,从根到所有物种的距离都相同。构建方法有很多,比较常用的就是r8s.

    这里用cafetutorial_prep_r8s.py构建r8s的批量运行脚本,然后提取超度量树

    python python_scripts/cafetutorial_prep_r8s.py -i twelve_spp_raxml_cat_tree_midpoint_rooted.txt -o r8s_ctl_file.txt -s 35157236 -p 'human,cat' -c '94'
    r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt
    tail -n 1 r8s_tmp.txt | cut -c 16- > twelve_spp_r8s_ultrametric.txt
    

    运行CAFE

    运行CAFE有两种模式,一种是CAFE的命令行模式,先执行cafe进行CAFE的shell, 然后在其中执行命令。另一种是脚本模式,也就是你先把命令编辑完成,然后用cafe script_to_run.sh运行。

    估计出生-死亡参数 \lambda

    CAFE的主要功能就是根据给定的进化树和基因家族数估计一个或多个 birth-death(\lambda)参数。\lambda 参数描述的是基因出现或者消失的概率。

    为整个树估计单个\lambda

    编辑cafetutorial_run1.sh。CAFE的命令不能有额外的空格出现在 tree后面的()中,以及lambda 的 -t 后的()中,否则运行时会无法正确解析文件导致报错。

    #!cafe
    load -i filtered_cafe_input.txt -t 4 -l reports/log_run1.txt
    tree ((((cat:68.710507,horse:68.710507):4.566782,cow:73.277289):20.722711,(((((chimp:4.444172,human:4.444172):6.682678,orang:11.126850):2.285855,gibbon:13.412706):7.211527,(macaque:4.567240,baboon:4.567240):16.056992):16.060702,marmoset:36.684935):57.315065):38.738021,(rat:36.302445,mouse:36.302445):96.435575)
    lambda -s -t ((((1,1)1,1)1,(((((1,1)1,1)1,1)1,(1,1)1)1,1)1)1,(1,1)1)
    report reports/report_run1
    

    然后运行如下命令

    mkdir -p reports
    cafe cafetutorial_run1.sh
    

    结果统计

    上一步运行结束后的报告文件在reports/reportrun1.cafe,可以用已有的脚本分析哪些基因家族发生了扩张或者搜索

    python  python_scripts/cafetutorial_report_analysis.py -i reports/report_run1.cafe -o reports/summary_run1
    

    reports文件夹下会出现如下文件

    • summary_run1_node.txt: 统计每个节点中扩张,收缩的基因家族数
    • summary_run1_fams.txt: 具体发生变化的基因家族

    为高基因拷贝数的基因家族预测\lambda

    之前过滤掉的高拷贝数变异的基因家族可以单独进行分析, 运行命令如下

    #!cafe
    load -i large_filtered_cafe_input.txt -t 4 -l reports/log_run2.txt
    tree ((((cat:68.710507,horse:68.710507):4.566782,cow:73.277289):20.722711,(((((chimp:4.444172,human:4.444172):6.682678,orang:11.126850):2.285855,gibbon:13.412706):7.211527,(macaque:4.567240,baboon:4.567240):16.056992):16.060702,marmoset:36.684935):57.315065):38.738021,(rat:36.302445,mouse:36.302445):96.435575)
    lambda -l 0.00265952 -t ((((1,1)1,1)1,(((((1,1)1,1)1,1)1,(1,1)1)1,1)1)1,(1,1)1)
    report reports/report_run2
    

    为树的不同部分预测多个\lambda

    如果你认为不同物种或者不同分支的基因家族进化速率不同,那么可以让CAFE预测多个\lambda值. 对lambda部分进行调整, 相同数字表示\lambda相同,不同数字表示\lambda不同。

    #!cafe
    load -i filtered_cafe_input.txt -t 4 -l reports/log_run3.txt
    tree ((((cat:68.710507,horse:68.710507):4.566782,cow:73.277289):20.722711,(((((chimp:4.444172,human:4.444172):6.682678,orang:11.126850):2.285855,gibbon:13.412706):7.211527,(macaque:4.567240,baboon:4.567240):16.056992):16.060702,marmoset:36.684935):57.315065):38.738021,(rat:36.302445,mouse:36.302445):96.435575)
    lambda -s -t ((((3,3)3,3)3,(((((1,1)1,2)2,2)2,(2,2)2)2,3)3)3,(3,3)3)
    report reports/report_run3
    

    CAFE主要输出文件格式

    CAFE主要输出内容如下,下游分析所需信息需要通过对其解析获得。

    2013053-12d6edceae9abf8d.png
    信息

    Lambda是整个进化树的预测值

    # IDs of nodes表示不同节点的编号,这里cat为0,horse为2,cat和horse所在的节点是1.

    最后是每个基因家族的结果。以最开始的表示行为例,第一列对应输入基因家族的编号;第二列是Newick的进化树,cat_59中的59表示该基因家族在cat里有59个基因;第三列是Family-wide P-value,用于表明该基因家族是否是显著性的扩张或是收缩,这里是0.438,说明变化不明显。在第三列的p值小于0.01时,第四列表明哪个分支的基因家族发生了变化,上图中只有ID 11的基因家族有变化, 但是0,1,2,4分支并没有变化。

    参考资料

    展开全文
  • RepeatModeler可用来从头对基因组的重复序列家族进行建模注释,它的核心组件是RECON和RepatScout。 使用方法 以拟南芥的参考基因组为例,假设基因组的名字为"Athaliana.fa" 第一步:为RepeatModeler创建索引数据库 ...

    RepeatModeler可用来从头对基因组的重复序列家族进行建模注释,它的核心组件是RECON和RepatScout。

    使用方法

    以拟南芥的参考基因组为例,假设基因组的名字为"Athaliana.fa"

    第一步:为RepeatModeler创建索引数据库

    BuildDatabase -name ath -engine ncbi Athaliana.fa
    # -engine ncbi: 表示使用rmblast
    # -name aht: 表示数据库的名字为ath
    

    这一步其实认为调用了makeblastdb,结果文章和makeblastdb一致

    第二步:运行RepeatModeler

    RepeatModeler -database ath -engine ncbi -pa 10 &> ath.out &
    # -database 要和上一步一致
    # -engine 要和上一步一致
    # -pa 表示线程数
    

    运行中的的文件存放在RM_.xxx文件夹下

    RM_82213.MonOct151054032018
    ├── consensi.fa 
    ├── consensi.fa.classified
    ├── consensi.fa.masked
    ├── families-classified.stk
    ├── families.stk
    ├── round-1
    ├── round-2
    ├── round-3
    ├── round-4
    └── round-5
    

    运行结束后,就得到了ath-families.faath-families.stk。 前者是找到的重复序列,后者是Stockholm格式的种子联配文件(seed alignment file), 可以用util/dfamConsensusTool.pl上传到Dfam_consensus数据库中。

    grep -i Unkown ath-families.fa
    # 没有结果,全都归类,毕竟拟南芥
    

    RM_.xxxconsensi.fa.classifiedath-families.fa内容一样,也是FASTA格式的文件,,只不过每个序列的ID会标注它来自于哪个重复序列家族,如果无法归类,则用"Unkown"标注。

    你可以将ath-families.fa分ModelerID.lib和Modelerunknown.lib,

    seqkit grep -nrp Unknown ath-families.fa > Modelerunknown.lib
    seqkit grep -vnrp Unknown ath-families.fa >  ModelerID.lib
    

    其中Modelerunknown.lib用RepeatMasker或TEclass进一步注释,如果能够被分类则从Modelerunknown.lib,移动到ModelerID.lib

    这里用TEclass进行分类,软件安装使用参考使用TEclass对TE一致性序列进行分类.

    TEclassTest Modelerunknown.lib
    

    这一步其实无法运行,因为没有输入。

    第三步:过滤基因片段

    许多文章都没有这样子做,因此这一步完全可以不用看

    RepeatModeler找到的重复序列进一步在植物蛋白数据库(不包括转座子蛋白)进行搜索,如果和植物蛋白匹配,或者在序列的侧翼50bp以内,就将该重复序列剔除。这一步可用工具是ProtExcluder
    ProtExcluder的安装需要HMMER3, 方法如下

    首先一定要装HMMER3.0,然后安装方法为

    wget http://eddylab.org/software/hmmer/hmmer-3.0.tar.gz
    tar xf hmmer-3.0.tar.gz
    cd hmmer-3.0
    ./configure --prefix=$HOME/opt/biosoft/hmmer-3.0
    make && make install
    cd easel
    ./configure --prefix=$HOME/opt/biosoft/hmmer-3.0
    make && make install
    

    后续安装的ProtExcluder1.1非常坑爹,他居然认为hmmer的运行文件是放在binaries目录下,所以你还需要去~/opt/biosoft/hmmer-3.0把bin文件夹改成binaries

    然后安装ProtExcluder1.1

    tar zxf  ProtExcluder1.1.tar.gz
    mv ProtExcluder1.1 ~/opt/biosoft
    cd  ~/opt/biosoft/ProtExcluder1.1 
    ./Installer.pl   -m   ~/opt/biosoft/hmmer-3.0 -p P ~/opt/biosoft/ProtExcluder1.1/
    # 注意运行路径后一定要有"/"
    

    然后将找到的重复序列用blastx比对到植物的蛋白数据库中,你可以到RefSeq上进行下载。

    blastx -query ath-families.fa -db /path/to/refseq-plant/db/plant.protein -num_threads 70 > ath-families.blastx &
    

    最后运行ProtExcluder.pl

    /path/to/ProtExcluder.pl -f 50   ath-families.blastx  ath-families.fa 
    

    从理论上说结果是"XXXnoProtFinal",但这个破软件各种报错,而且报错信息信息量太少,直接放弃这个破软件了。

    考虑到目前看了很多文献也没人说要过滤,所以这一步就放弃吧

    注1: 一般而言RepeatModeler在一周内就能运行完,物种小一点,服务器好一点,基本上一天就行了
    注2: 一般而言,同源注释所用的重复序列库的数据比较可靠,但不一定全,而我们自己建立的重复序列库比较全,但是未必准。
    注3: SwissProt的植物蛋白数据库的下载地址为:ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions,而NCBI的Refseq植物蛋白数据库下载地址为ftp://ftp.ncbi.nih.gov/refseq/release/plant/


    版权声明:本博客所有文章除特别声明外,均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。

    2013053-24d1e185d725662e.png
    扫码即刻交流
    展开全文
  • 提取基因序列文件

    2018-11-28 14:51:32
    perl可编译代码,可用于提取目的序列文件中所需序列id
  • 如何对基因组序列进行注释

    万次阅读 多人点赞 2018-09-07 18:01:37
    基因组组装完成后,或者是完成了草图,就不可避免遇到一个问题,需要对基因组序列进行注释。注释之前首先得构建基因模型,有三种策略: 从头注释(de novo prediction):通过已有的概率模型来预测基因结构,在预测...
  • 写在前面通过我几天的学习,我发现,seqkit十分好用,将序列的各种操作都囊括进去,加入多线程,我个人认为这将是非常好的胶水,在处理无论是基因组还是其他组学。定是一个必学神器。注意一下教程...
  • WGDI(全基因组重复识别),一种基于 Python 的命令行工具,可让研究人员深入了解递归多倍化并进行跨物种基因组比对分析。 官方文档 下一 篇会选择一个物种的分析结果做示例。 安装 ## 1.使用conda安装 conda ...
  • 基因组共线性工具MCScanX使用说明

    千次阅读 2018-11-21 15:09:31
    基因组共线性工具MCScanX使用说明 简介 MCScanX工具集对MCScan算法进行了调整,用于检测共线性和同线性区域,还增加了可视化和下游分析。。MCscanX有三个核心工具,以及12个下游分析工具 实际使用 第一步: BLASTP比对...
  • 我梳理了GWAS全基因组关联分析的整个流程,并提供了基本的命令,用到的软件包括BWA、samtools、gatk、Plink、Admixture、Tassel等,在此分享出来给大家提供参考。 一、BWA比对 1.构建索引 bwa index -a is ...
  • 我梳理了GWAS全基因组关联分析的整个流程,并提供了基本的命令,用到的软件包括BWA、samtools、gatk、Plink、Admixture、Tassel等,在此分享出来给大家提供参考。 一、BWA比对 1.构建索引 bwa index -a is example....
  • awk 'BEGIN{OFS=FS="\t"}ARGIND==1{id[$1]=1}ARGIND==2{if(id[$1]==1) print $0}' GRCh38_tmp.seqid GRCh38.gtf >GRCh38_tmp.gtf # 统计构建基因组索引所需的计算资源 /usr/bin/time -o STAR_human_genome_${i}.log ...
  • 先将gff转成bed格式, python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_thaliana.TAIR10.44.gff3.gz > ath.bed python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id ...
  • $ seqkit stat input.fasta file format type num_seqs sum_len min_len avg_len max_len input.fasta FASTA DNA 1 293,737 293,737 293,737 293,737 用 minimap2 将原来的序列和组装结果进行比对,并按照 「R...
  • 第一次进行基因数据的处理 将测序出来的数据构建一个prifile 关于节点篇 一般的服务器集群都会分为登录节点和计算节点,为了能够选择到计算节点首先了解节点的情况 使用SSH切换到计算节点 ssh cngb-9 #cngb-9为一个...
  • 2、根据基因ID从fasta文件中提取序列序列 3、比对后的fasta格式转化为phylip格式 4、提取四倍兼并位点 5、根据基因名字从gtf或者gff文件中提取基因ID 二、可以模仿一些常用的软件的部分功能,去写能够实现相同功能的...
  • 已知,我们通过seqkit faidx ref.fa chr8:25234310-25266151 > target.fa提取基因组上的一个片段序列,接着我们用AUGUSTUS对这段序列进行预测 augustus --species=arabidopsis target.fa > target.gff 我们的...
  • 如果一个物种在演化过程中发生过多倍化,那么在基因组上就会存在一些共线性区域(即两个区域之间的基因是旁系同源基因,其基因的排布顺序基本一致)。 例如拟南芥经历了3次古多倍化,包括2次二倍化,一次3倍化 (Tang. ...
  • 如何用WGDI进行共线性分析(上)

    千次阅读 2021-02-01 14:50:09
    如果一个物种在演化过程中发生过多倍化,那么在基因组上就会存在一些共线性区域(即两个区域之间的基因是旁系同源基因,其基因的排布顺序基本一致)。 例如拟南芥经历了3次古多倍化,包括2次二倍化,一次3倍化 (Tang. ...
  • 前面我们评估了不同大小基因组基于STAR构建索引所需的计算资源和时间资源、不同大小数据集基于STAR进行比对所需的计算资源和时间资源和STAR比对速度与分配线程的关系。将人类基因组按染色体...
  • 我首先是看基因家族鉴定的文献,发现很多是用 domain建树。如果此时还是不好,那么树是什么样,就是什么样。 首先要获取对应的转录家族序列,这一步简单,只要去pfamA看一下,但是一定要尽量细化到最下面的亚家族。...
  • 如何从BAM文件中提取fastq

    千次阅读 2018-02-12 23:01:00
    虽然高通量测序分析最常用的操作是将fastq比对到参考基因组得到BAM文件,但偶尔我们也需要提取BAM文件中特定区域中fastq。最开始我认为这是一个非常简单的操作,因为samtools其实已经提供了相应的工具samtools fastq...
  • 基因组扩增子最新分析流程QIIME2-了解分析趋势 ,读过的朋友会发现,里面的每个分析流程中都不再使用聚类方法生成OTU,而是调用DADA2 [1]对原始数据进行去噪,相当于以100%的相似度聚类,而仅仅对低质量序列进行...

空空如也

空空如也

1 2 3
收藏数 43
精华内容 17
关键字:

seqkit根据基因id