cellranger参考数据库gtf过滤（替换和删除）

cellranger在建库之前，有时候需要对GTF文件进行编辑修改。

比如为了后续在Seurat中使用正则匹配去除线粒体的基因，就需要修改原来GTF中的线粒体基因symbol名称，比如加上MT-或mt-前缀。

比如发现有些基因组存在一些重复的ID或者symbol的情况：

GCF_000003025.6_Sscrofa11.1_genomic.gtf筛选部分结果

这种情况需要重命名其中一个基因symbol，也就是需要在GTF里面进行替换。

比如还有一些基因，相同的基因ID，但是却在不同的染色体上：

GCF_000003025.6_Sscrofa11.1_genomic.gtf筛选部分结果，分别表示染色体和基因ID

这种情况需要删除掉其中一条描述信息，也就是需要在GTF里面进行删除。

上面几种情况，要根据GTF里面具体的第九列中的描述内容针对性地进行替换。
先放出来处理的Shell脚本如下：

#!/bin/bash
#!/usr/bin/env bash 
set -e
# echo the help if not input all the options
help()
{
cat <<HELP
---------------------------------------------------------------
     Author: Myshu                                                                            
     Version: 1.0                                              
     Date: 2022/12/29
     Description: GTF文件中替换gene symbol名称或删除重复ID的不同chr的行
---------------------------------------------------------------
USAGE: $0 <dup File> <GTF File> <OUT File> <Tag>
    or $0 -h # show this message
EXAMPLE:
    $0
NOTE
  1.Dup File: 重复ID的文件，两列。分别表示基因ID + symbol或者ID + chr 
  2.GTF File：需要替换的GTF文件（要注意关键标注是否一致）
  3.默认当前路径下必须有校正后的基因ID + symbol + 所在的染色体三列的对应表reference.xls
  4.Tag: 可选 mit/symbol/chr，分别表示替换线粒体基因名称，重名的gene symbol和删除重复ID的chr行
HELP
exit 0
}
[ -z "$1" ] && help
[ "$1" = "-h" ] && help

file=$1
gtf=$2
out=$3
tag=$4

echo ========== Start at : `date "+%Y-%m-%d/%H:%M:%S"` ==========

IFS_old=$IFS
IFS=$'\n'
> sed.$out
for l in `cat $file`
do
    id=`echo $l |cut -f 1 | tr '\n' ' ' | sed -e 's/[ ]*$//g'`
    name=`echo $l |cut -f 2 | tr '\n' ' ' | sed -e 's/[ ]*$//g'`
    #echo $id $name
    if [ $tag == "symbol" ];then
    res=`grep -P "^$id\t" reference.xls | cut -f 2 | tr '\n' ' ' | sed -e 's/[ ]*$//g'`
    #echo $res  # gene symbol name
    if [ $name != $res ];then
      echo -e "$id\t$name\t$res"
      #sed -i "s/\(.*gene_id \"$id\";.*gene_name \"\)[^\"]*\(\";.*\)/\1$res\2/g" $out
      echo -e "s/\(.*gene_id \"$id\";.*gene_name \"\)[^\"]*\(\";.*\)/\1$res\2/g" >> sed.$out
    fi
  elif [ $tag == "mit" ];then
    res=`grep -P "^$id\t" reference.xls | cut -f 2 | tr '\n' ' ' | sed -e 's/[ ]*$//g'`
    #echo $res
    if [ $name != $res ];then
      echo -e "$id\t$name\t$res"
      #sed -i "s/\"$name\"/\"$res\"/g" $out
      echo -e "s/\"$name\"/\"$res\"/g" >> sed.$out
    fi
  else
    chr=$name
    res=`grep -P "^$id\t" reference.xls | cut -f 3 | tr '\n' ' ' | sed -e 's/[ ]*$//g'`
    #echo $res  # chr name
    if [ $chr != $res ];then
            echo -e "$chr\t$id\t$res"
            #sed -i "s/^$chr.*gene_id \"$id\".*$//g" $out
            echo -e "s/^$chr.*gene_id \"$id\".*$//g" >> sed.$out
    fi
  fi
done;
IFS=$IFS_old
sed -f sed.$out $gtf > $out
echo ========== End at : `date "+%Y-%m-%d/%H:%M:%S"` ==========

脚本中主要针对上面描述的几种情况进行了替换和删除。

• 要注意这个脚本并不是对所有GTF都通用，我们在实际分析中，需要确定好实际的GTF第九列中的关键词和格式，再进行编辑脚本中的sed命令行代码，才能适用。
• sed命令跑的有点慢，建议放后台，我一个基因组2-3万多个基因大概花了8个多小时。大家如果有可以加速的也可以给出建议_我再改进
• 另外，别忘了替换之后的校验，避免遗漏了一些GTF中的行或者替换错误的情况发生！

欢迎大家批评指正！