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中的行或者替换错误的情况发生！
欢迎大家批评指正！
网友评论