对某个物种进行分子进化分析，保守序列提取以及CRISPR-gRNA设计

基因组序列的获取

这里我们直接获取完整的基因组组装序列
在ncbi的基因组数据中检索物种关键词，获取该物种所有组装好的序列，去FTP下载

• 基因组fasta序列（*.fna)
• 基因组注释文件(*.gff)

根据gff文件获取每个想要基因的序列

基因的获取：
通过查找相关物种的文献，提取出该物种的有效基因
从gff文件中获取位置坐标并提取序列
原理很简单，就是按照小标题所示，可以自己写一个脚本实现，这里我用python实现

import argparse
import os
import glob
import re
import pandas as pd
from Bio import SeqIO
from Bio.Alphabet import IUPAC
from Bio.Seq import Seq
parser = argparse.ArgumentParser(
    description='This script was used to get seq from gff file')
parser.add_argument('-g', '--gtf', help='Please input your bam file')
parser.add_argument('-i', '--id', help='Please input your reference file')
#parser.add_argument('-f', '--fasta', help='Please input start position')
args = parser.parse_args()


def read_gtf(file):
    with open(file) as f:
        alldata = f.readlines()
        start = pd.Series([line.split("\t")[3]
                           for line in alldata if not line.startswith("#") and line.split("\t")[2] == "gene"])
        stop = pd.Series([line.split("\t")[4]
                          for line in alldata if not line.startswith("#") and line.split("\t")[2] == "gene"])
        strand = pd.Series([line.split("\t")[6]
                            for line in alldata if not line.startswith("#") and line.split("\t")[2] == "gene"])
        genename = pd.Series([re.findall("Name=\w+", line.split("\t")
                                         [8])[0].replace("Name=", "") for line in alldata if not line.startswith("#") and line.split("\t")[2] == "gene"])
    gtf_data = pd.DataFrame(
        {'start': start, 'stop': stop, 'strand': strand, 'genename': genename})
    # print(gtf_data.head)
    return gtf_data


def get_gtffile(path):
    name = glob.glob(os.path.join(path, "*.gff"))
    # print(name)
    return name


def read_fasta(file):
    for seq_record in SeqIO.parse(file, "fasta"):
        fasta = str(seq_record.seq)
    return fasta


def read_fasta2(file):
    for seq_record in SeqIO.parse(file, "fasta"):
        nameid = seq_record.id.split(" ")[0]
    return nameid


def main():
    sequence = {}
    gtf_file_lst = get_gtffile(args.gtf)
    id = args.id
    i = 0
    for gtf in gtf_file_lst:
        gtf_data = read_gtf(gtf)
        if not gtf_data.empty:
            all_sequence = read_fasta(".".join(gtf.split(".")[0:-1])+".fna")
            all_sequence_name = read_fasta2(
                ".".join(gtf.split(".")[0:-1])+".fna")
            if not gtf_data[gtf_data['genename'] == id].reset_index().empty:
                start = gtf_data[gtf_data['genename']
                                 == id].reset_index().loc[0, 'start']
                stop = gtf_data[gtf_data['genename']
                                == id].reset_index().loc[0, 'stop']
                strand = gtf_data[gtf_data['genename'] ==
                                  id].reset_index().loc[0, 'strand']
                # print(start)
                # print(stop)
                # print(strand)
                transtable = str.maketrans("AGCT", "TCGA")
                if strand == "+":
                    sequence[">{}_{}".format(all_sequence_name,
                                             id)] = all_sequence[int(start)-1:int(stop)]
                else:
                    sequence[">{}_{}".format(all_sequence_name,
                                             id)] = all_sequence[int(start)-1:int(stop)][::-1].translate(transtable)
            i = i+1
    with open("{}.fasta".format(id), 'a') as f:
        for key, val in sequence.items():
            f.write(key+"\n"+val+"\n")


if __name__ == "__main__":
    main()

直接输入-g gtf文件夹位置 -i genename

对所有的基因序列做保守序列分析

我这里使用了meme软件自动发现高度保守模型
可以使用在线版本，但并不适合大量序列
所以直接使用离线版

meme coding/*.fasta -dna -oc . -mod zoops -nmotifs 3 -minw 20 -maxw 100 -objfun classic -minsites 2 -markov_order 0

因为后续要设计gRNA,
所以限制序列长度在20-100bp之间（-minw,-maxw)，
zoops模式是meme发现模式中保守型和速度相对居中的一种模型(-mod)，
-oc代表新生成的文件会覆盖原来的

接着写脚本提取保守序列

import argparse
parser = argparse.ArgumentParser()
parser.add_argument("-i", "--input", help="increase output verbosity")
args = parser.parse_args()
with open("meme.txt") as f:
    i=1
    for line in f:
        if line.startswith("MOTIF"):
            print(">{}_{}".format(args.input,i))
            print(line.split(" ")[1])
            i=i+1

gRNA设计

大规模gRNA设计建议采用FlashFry
建立基因组索引文件

mkdir tmp
java -Xmx4g -jar FlashFry-assembly-1.9.2.jar \
 index \
 --tmpLocation ./tmp \
 --database Hg19_database \
 --reference Hg19.fa \
 --enzyme spcas9ngg

寻找gRNA序列

java -Xmx4g -jar FlashFry-assembly-1.9.2.jar \
 discover \
 --database hg19_database \
 --fasta A179L.motif.fasta \
 --maxMismatch 3 \ #建议设置该参数
 --output A179L.motif.output

分析脱靶个数并进行排序，提取出想要的序列

import pandas as pd
import glob

#我这里分别对人类基因组和猪基因组进行了gRNA设计并合并otcount,排序
def get_gRNA(gene):
    df1=pd.read_csv("{}.motif.output".format(gene),header=0,sep="\t")
    df2=pd.read_csv("{}.motif.output2".format(gene),header=0,sep="\t")
    df1=df1[['contig','target','otCount']].set_index('target')
    df1.columns=['contig','otCount_human']
    df2=df2[['contig','target','otCount']].set_index('target')
    df2.columns=['contig','otCount_pig']
    fin=df1.join(df2[['otCount_pig']])
    fin=fin.sort_values(['otCount_pig','otCount_human']).iloc[:3,]
    return fin

old_table=get_gRNA("A179L")
#这里随便找个基因作为初始化DataFrame
for gene in glob.glob("*.output"):
    new_gene=gene.replace(".motif.output","")
    #print(new_gene)
    table=get_gRNA(new_gene)
    #print(table)
    old_table=pd.concat([table,old_table])
#print(old_table.head)

old_table.to_csv("coding2.gRNA.csv",sep="\t")
~            

构建物种树

为了结合多个基因的结果，并且方便建树，采用了ETE3

基于几个连锁比对重建物种树需要选择两个工作流程：
i）用于比对每个基因家族序列的基因树工作流程（-w）
ii）基于超矩阵连接和构建树的工作流程alignment（-m）

所有基因/蛋白质的序列必须在单个fasta文件中传递，并且还需要COG（直系同源组）文件。COGs文件必须是包含与输入文件中相同的序列ID的文本文件。每个TAB分隔线将被视为COG。
例如，以下示例将分别定义3个大小为3,2和4的COG：

sp1_seqA sp2_seqA sp3_seqA
sp1_seqB sp2_seqB    
sp1_seqC sp3_seqC sp4_seqC sp5_seqC

建树命令：

ete3 build -a proteome_seqs.fa --cogs cogs.txt -o sptree1_results -m sptree_fasttree_100 -w standard_fasttree

这里在安装ETE3的时候千万不要使用python3,有极大的可能会莫名其妙报错
如果想要覆盖已有的输出，使用--clearall
cogs文件实例，蛋白序列实例

所以如果想要正确使用的话，我建议使用miniconda创建python2的虚拟环境进行安装

conda create -n ETE python=2
source activate ETE
conda install -c etetoolkit ete3 ete_toolchain

NUP62.aa.fa.final_tree.png