靶标序列设计

具体参数如下

-i <str>                 Input file
-x <int>                 Length of sgRNA[20]
-l <int>                 The minimum value of GC content [20]
-m <int>                 The maximum value of GC content [80]
-g <str>                 The reference genome sequence
-o <str>                 Searching CRISPR target sites using DNA strands based option(s/a/b)
                         [s, sense strand searching mode]
                         [a, anti-sense strand searching mode]
                         [b, both strand searching mode]
-t <str>                 Type of sgRNA searching mode(s/p)
                           [s, single-gRNA searching mode]
                           [p, paired-gRNA searching mode]
-v <str>                 Operation system(w/l/u/m/a)
                            [w, for windows-32, 64]
                            [l, for linux-64]
                            [u, for linux-32]
                            [m, for MacOSX-64]
                            [a, for MacOSX-32]
-n <int>                 Maximum number of mismatches [5]
-s <int>                 The minimum value of sgRNA offset [-2]
-e <int>                 The maximum value of sgRNA offset [32]
-p <str>                 Output path 

一般性的命令

需要注意的是，在输入的CDNA序列中gene的ID编号不能够带有小数点，并且这里p参数后面可以省略参数，不然老是报错

perl sgRNAcas9_3.0.5.pl -i genes_end.fasta -x 20 -l 40 -m 60 -g Ghirsutum_genome.fasta -o b -t s -v l -n 5 -s -3 -e 33 -p /public/home/llpei/zpliu/zpliu/software/sgRNAcas9_3.0.5 2>log.txt

将所有OT文件整合在一个文件中

cat A.Sort_OT_byID/*  >all_genen_OT.txt

提取对应等级的靶标ID

Discard > High_risk > moderate_risk > low_risk > repeat_sites_or_bad ? > Best
一般提取Best 、repeat_sites_or_bad ? 、low_risk等级的ID

0M    1M      2M    3M     rank 
2       0       0     0    repeat_sites_or_bad
1       0       0    5    low_risk
1       0       0    3    Best 

例如repeat等级中0M靶向的位置有2个，我们要看看它靶向的位置是否是同一个基因，进行sgRNA评价

//提取Best靶标中0个错配的sgRNA 信息
awk -F "\t" '$23~/^B/{print $1}' sgRNAcas9_report.xls  |xargs  -I {} grep {} all_genen_OT.txt |awk -F "\t" '$4~/0M/{print}' >Best_OT.text
//提取repeat_sites_or_bad? 类型的信息
awk -F "\t" '$23~/^R/{print $1}' sgRNAcas9_report.xls  |xargs  -I {} grep {} all_genen_OT.txt |awk -F "\t" '$4~/0M/{print}' >repeat_OT.text
//提取低质量low_risk类型的信息
awk -F "\t" '$23~/^L/{print $1}' sgRNAcas9_report.xls  |xargs  -I {} grep {} all_genen_OT.txt |awk -F "\t" '$4~/0M/{print}' >Low_OT.text

靶标序列评价

使用自带的脚本ot2gtf_v2.pl，对得到的sgRNA靶标进行评价，看这个靶标是否靶向哪个基因的外显子区域,

在OT输出文件中包含了sgRNA在基因组上的位置信息，所以想借助gtf文件来看一下对应的位置是不是基因的exon区域，并且将对应的基因名字添加到后面一栏，基因组的gtf文件可以借助其他脚本进行转换获取

Ghir_A05G009820_A_54    POT984232       - - - - - - - - - - - - - - - - - - - - - - -   0M      Ghir_A05        8868917 +
Ghir_D05G008800_S_140   POT981667       - - - - - - - - - - - - - - - - - - - - - - -   0M      Ghir_A05        8060310 -
Ghir_D05G008800_S_140   POT3810799      - - - - - - - - - - - - - - - - - - - - - - -   0M      Ghir_D05        7188361 -
Ghir_D05G008800_S_141   POT981663       - - - - - - - - - - - - - - - - - - - - - - -   0M      Ghir_A05        8060259 -

• 执行命令

  perl ../Usefull_Script/ot2gtf_v2.pl  -i Low_OT.text  -g ../Gh_gene.gtf  -o Low_OT_gtf_out.text
  也可以将上一步所有的OT文件合并之后，在运行脚本
  

这个脚本50行左右开始需要修改一下，默认是使用人类染色体编号进行的

    #1    pseudogene    gene    11869    14412    .    +    .    gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene";

    (my $gtf_chr, my $gtf_annotation1, my $gtf_annotation2,    my $gtf_beg, my $gtf_end, 
        my $get_dot1, my $gtf_strand, my $gtf_dot2, my $gtf_annotation3)=split/\t/, $_; 

    # next unless $gtf_chr=~ /^[0-9]|^[X]|^[Y]/;

    next unless $gtf_chr eq $chr;

    # next if !$gtf_beg =~ /\d+/;
    #gene_id "ENSG00000223972"; transcript_id "ENST00000515242"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-201"; transcript_source "ensembl";
    if ($gtf_annotation3=~ /".*gene_name "([^"]+)"/) {

    $geneID = $1;
    #print $geneID,"\n";

去除脱靶的sgRNA

从三个点出发

1. 靶标基因与sgRNA的序列信息一致赋权值 0
2. sgRNA靶向同源基因和自己本身赋权值 0
3. 靶标序列脱靶赋权值 1
4. 最后将同一个sgRNA靶标的权值相加，为0则表示没有脱靶；否则脱靶舍弃

#awk脚本
-F "\t" '{
    a=substr($1,1,15);a1=substr($1,6,1);a2=substr($1,7,2);b1=substr($5,6,1);b2=substr($5,7,2)
    }{
if(a==$8)print $0"\t"0;
else if(a1=="A"&&b1=="D"){
    if(a2==b2||a2=="02"&&b2==03||a2=="03"&&b2=="02"){
        print $0"\t"0;}else print $0"\t"1;}
        else if(a1=="D"&&b1=="A"){
            if(a2==b2||a2=="02"&&b2=="03"||a2=="03"&&b2=="02"){
                print $0"\t"0;}else print $0"\t"1;}
else print $0"\t"1;}' Best_Repeat_Low_OT_gtf  >2222222

从上一步结果中提取没有脱靶的sgRNA信息

awk -F "\t" '{a[$1]+=$9}END{for(i in a){if(a[i]==0)print i}}' 2222222 |xargs  -I {} grep -e "{}\s" Best_Repeat_Low_OT_gtf  >sgRNA_end

得到我们所需的sgRNA靶标序列

对合格靶标序列进行筛选

经过上一步筛选后的sgRNA文件，我们需要根据以下几个指标筛选比较理想的靶标序列

1. 靶标序列尽量靠近5’端
2. 同一个基因找两个靶标序列，尽量让这两段序列间隔在100bp左右

写了达到目的写了一个python脚本，脚本的使用如下

usage:

    -h|--help    print help information
    -g|--gff=    gff file path way
    -s|--sgRNA=    sgRNA file path way
    -l|--genelength=    length of gene
    -r|--sequence=    sgRNA sequence path way
    -o|--outfile=    output file path way

脚本获取

参考

tiramisutes
sgRNAcase9awk