Без опису

hdb cd53d00105 test		8 місяців тому
script	cd53d00105 test	8 місяців тому
Readme.md	ca14764da2 first commit	8 місяців тому

真核有参参考基因组整理

参考基因组目录命名规范：
只包含字母数字. _
有数据库：拉丁名_数据库名称_版本例：mus_musculus_ensembl_111 Homo_sapiens_NCBI_GCF_000001405.40_GRCh38.p14
无数据库（客户提供）：拉丁名_合同号_客户姓名全拼：例： Pseudomonas_fluorescens_PN20240113005_liuyang

非常规gff纠错和矫正为ensembl格式gtf

先修复gff为标准的gff3格式；然后使用genpred转为gtf

gff 转为 gtf

gff3ToGenePred $gff genpred
genePredToGtf file genpred example.gtf
# 下面这一步会自动加上gene_biotype "protein_coding"；请根据gff文件里是否有基因类型的信息；去处理
perl /Business/psn_company/t04/hudabang/bin/yaogu/tools/add_gene_type_to_gtf.pl -i example.gtf -o  example.add_gene.gtf
# 注意这个里面是否有 gene_biotype 信息；没有的话需要对有cds的基因加上 gene_biotype "protein_coding"
sed '/\tgene\t/s/$/ gene_biotype "protein_coding";/' example.add_gene.gtf

方法一：使用工具AGAT

source /PERSONALBIO/work/transcriptome/t00/software/biosoftware/miniconda3/bin/activate AGAT_gff
export PERL5LIB=/PERSONALBIO/work/transcriptome/t00/software/biosoftware/miniconda3/envs/pfam_scan/share/pfam_scan-1.6-4
# 修复gff 缺少基因或mRNA;cds没有对应的外显子 
agat_convert_sp_gxf2gxf.pl  -g example.gff -o example.fixed.gff3
# gff3 转gtf ；参考上面

方法二：手动处理；具体问题具体分析

常见问题：
gene cds exon mRNA 类型不全
CDS feature 第8列不为 0 1 2
Parent和ID相同
正负链为. 或者同一个基因下出现不同的 + -
出现不同位置；相同的基因id
。。。待补充

参考基因组整理自动化流程

/Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/main.sh -h

Usage: /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/main.sh [-h|--help] [-i|--indir <indir>] [-o|--outdir <outdir>] [-c|--config <config>] [-n|--cpu <cpu>] [-M|--memory <memory>] [-v|--version] [-a|--allgene]
Options:
  -h, --help            Show this help message and exit.
  -i, --indir           Input directory containing the files that in config file
  -o, --outdir          Output directory
  -c, --config          Config file for the pipeline.
  -n, --cpu             Number of CPUs in bsub;default: 6
  -M, --memory          Memory limit  in bsub;default: 50G
  -v, --version         Show the version of the pipeline and exit.
  -a, --allgene         Include all genes, not just protein-coding genes.

#示例路径
/Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/example/mus_musculus_ensembl_111
#示例命令
/Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/main.sh -i /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/example/mus_musculus_ensembl_111 -o /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/example/mus_musculus_ensembl_111/genome -c  /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/example/mus_musculus_ensembl_111/genome_conf.ini

参数说明

-i 输入路径；下面包含需要整理基因组用到的文件；基因组序列和基因结构注释gtf；蛋白序列或cds序列（序列title必须是gtf里的基因id）；其他功能注释文件等
-o 输出目录；
-c 基因组配置文件
-n 分析过程中比对任务投递到集群的 -n
-M 分析过程中比对任务投递到集群的 -M ；默认50G；因为eggnogmap 至少需要50G
-v 流程版本
-a 整理所有基因(医学需要所有基因)还是protein_coding（真核有参mRNA）;默认 protein_coding

配置文件 genome_conf.ini 介绍

#详细解释见 ：/Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/genome_conf.ini
latin_name=Mus_musculus                                             物种拉丁名
taxid=10090                                                         物种taxid
url=https://ftp.ensembl.org/pub/release-111/fasta/mus_musculus/dna/ 物种数据库下载链接
version=111                                                         数据库上对应物种版本
dna_fasta=Mus_musculus.GRCm39.dna.toplevel.fa                       基因组序列
gtf=Mus_musculus.GRCm39.111.gtf                                     基因组结构注释gtf
fasta=pep.fa                                                        蛋白序列或cds序列（序列title必须是gtf里的基因id）
fasta_type=proteins                                                 序列类型：参数 proteins 或 CDS
id_table=id_table.txt                                               gtf里的基因id locus_tag old_locus_tag entry_id 蛋白id 基因类型
kindom=animal                                                       物种大类：参数： animal plant fungi
org=mmu                                                             物种kegg三字符； 没有则写物种大类
go_file=GO.txt                                                      GO注释文件
kegg_file=none                                                      KEGG注释文件
eggnog_file=none                                                    eggnog注释文件
swissprot_file=Swissprot.txt                                        swissprot注释文件
nr_file=NR.txt                                                      NR注释文件
uniprot_tsv=none                                                    uniprot 本物种的 GO注释
uniprot_fasta=none                                                  uniprot 本物种的 蛋白序列

参考基因组整理分步流程

参考基因组需要整理的文件 refCfg.ini 说明

DNA=DNA/Mus_musculus.GRCm39.dna.primary_assembly.fa 基因组序列和索引
RefDatabase_info=DNA/RefDatabase_info 基因组信息
chr_list=DNA/chr_list 染色体
hrds=DNA/hrds.txt hrds
cytoBand=DNA/cytoband.txt 圈图使用的染色体信息
GTF=GTF/Mus_musculus.GRCm39.108.gtf 结构注释gtf
GTFbed=GTF/GTF.bed 结构注释bed格式
hisat2SS=GTF/genome.ss 比对和可变剪切使用
lncRNA_gtf=GTF/lncRNA.gtf lncRNA结构注释gtf
DEXSeq_gff=GTF/DEXSeq.gff 差异外显子使用的外显子结构注释gff
matureMiRNA=mmu_mature.fa 已知成熟miRNA序列
annoDB=annoDB SNP annova使用
ppi=PPI/PPI.txt 蛋白互作文件
Annotation=Annotation/Annotation.xls 基因功能注释
Anno_summary=Annotation/Annotation.stat 基因功能注释统计
09_TFs 转录因子；真菌无转录因子

DNA相关文件

# 2.1 RefDatabase_info文件：
Genome	Homo_sapiens.GRCh38.dna.primary_assembly.fa
Path	http://asia.ensembl.org/Homo_sapiens/Info/Index
Genebuild by	Ensembl
Database version	104.38
Base Pairs	3,099,750,718

统计Base Pairs可用以下命令：
perl /Business/psn_company/t01/public/ScriptManager/lib/PrepRef/calculate_fasta_basepair/calculate_fasta_basepair.pl -i reference.fa

# 2.2 对DNA建立 hisat2 索引和fai索引
hisat2-build -p 12 reference.fa reference.fa
samtools faidx reference.fa

# 2.3 chr_list文件
chr_list 主要是用于画基因组圈图的，指定哪些染色体需要画到圈图中，一般为常染色体加性染色体，不包括没拼好的 scalfold，每行写一个染色体；若全是contig或scaffold, 则全部保留；从参考基因组序列文件中获取。
grep "^>" reference.fa | grep chromosome | cut -d " " -f1 | sed 's/^>//' > chr_list

# 2.4 生成 hrds 文件（SNP分析中用到）
该文件记录了物种和每个染色体长度。
# 参数 -org 是给物种拉丁名，下面以人为例：
perl /Business/psn_company/t01/public/ScriptManager/lib/PrepRef/create_hrds/create_hrds.pl -i reference.fa -o hrds.txt -org Homo_sapiens

# 2.5 cytoBand文件（用于基因组圈图绘制）
perl /Business/psn_company/t01/public/ScriptManager/lib/PrepRef/create_cytoband/create_cytoband.pl -i hrds.txt -c chr_list -o cytoband.txt

GTF相关文件

# 3.1 生成genome.ss文件 (hisat2 比对中用于寻找剪切位点的文件):
extract_splice_sites.py reference.gtf > genome.ss

# 3.2 生成 GTF.bed 文件 (bed 文件记录了基因的位置信息，主要用于 mapQC 的几个模块)
/Business/psn_company/t01/public/ScriptManager/lib/PrepRef/gtf2bed/gtf2bed.pl reference.gtf > GTF.bed

# 3.3 DEXSeq_gff (用于差异外显子分析)
python /Business/psn_company/t01/public/ScriptManager/lib/DE/exon_DEXSeq/dexseq_prepare_annotation.py reference.gtf DEXSeq.gff

# 3.4 lncRNA.gtf (此文件是 lncRNA 的结构注释文件，从GTF中提取，如果GTF中没有，下述命令会生成一个空文件)
perl /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/abstract_lncRNA_gtf.pl -i reference.gtf -o lncRNA.gtf

annoDB文件夹

mkdir annoDB
gtfToGenePred -genePredExt reference.gtf annoDB/Anno_refGene.txt
perl /Business/psn_company/t01/public/software/annovar/retrieve_seq_from_fasta.pl --format refGene --seqfile reference.fa annoDB/Anno_refGene.txt --out annoDB/Anno_refGeneMrna.fa

基因功能注释

GO注释

# 优先级1：NCBI 和ensembl 等公共数据库上 下载 
# 优先级2：uniprot上下载 对应和比对
wget -c -t 3 -q https://ftp.uniprot.org/pub/databases/uniprot/knowledgebase/reference_proteomes/Eukaryota/UP000005640/UP000005640_9606.fasta.gz
wget -c -t 3 -q https://ftp.uniprot.org/pub/databases/uniprot/knowledgebase/reference_proteomes/Eukaryota/UP000005640/UP000005640_9606.dat.gz
perl /personalbio1/R_Report/Users/hudb/script/download_genome/uniprot_dat2tsv.pl UP000005640_9606.dat.gz > UP000005640_9606.tsv
## 对应：
UP000005640_9606.tsv里分别是Uniprot AC ;GeneID (NCBI Entryid); GO ID 
## 比对
/Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/blast_uniprot.sh -i pep.fa -d UP000005640_9606.fasta -g UP000005640_9606.tsv
Usage:
        /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/blast_uniprot.sh -i cds.fa -o outdir -t CDS -d uniprot.fa  -g uniprot.go.tsv
OPTIONS:
        -h: help information
        -i: fasta file
        -o: output dir
        -t: proteins or CDS default:proteins
        -d: uniprot species fasta
        -g: uniprot species go annotation

# 优先级3：eggnog_map 注释
/Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/eggNOG_GO.sh -i pep.fa -s animal -o ./ -t proteins

Usage:
        /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/eggNOG_GO.sh -i cds.fa -o outdir -t CDS -s animal
OPTIONS:
        -h: help information
        -i: fasta file
        -o: output dir
        -t: proteins or CDS default:proteins
        -s: species [ animal plant fungi prokaryotes ]

KEGG注释

# 优先级1：脚本内会判断是否匹配上已知数据库
Usage:
        /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/KEGG.sh -i cds.fa -o outdir -t CDS -s animal -l idlist
OPTIONS:
        -h: help information
        -i: fasta file
        -o: output dir
        -l: idlist;col1-6 is gene_id  locus_tag    old_locus_tag  entry_id protein_id  gene_name
        -t: proteins or CDS default:proteins
        -s: kegg org [ animal plant fungi prokaryotes hsa mmu ...  ]
        -r: rate of Matched to known KO; default:0.5

eggnog注释

# 优先级1：NCBI 和ensembl 等公共数据库上 下载 
# 优先级2: eggnog_map 注释
同GO注释优先级3；会同时出eggnog 和GO注释结果

Swissprot注释

# 优先级1：NCBI 和ensembl 等公共数据库上 下载 
# 优先级2：比对
Usage:
        /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/Swissprot.sh -i cds.fa -o outdir -t CDS
OPTIONS:
        -h: help information
        -i: fasta file
        -o: output dir;default:./
        -t: proteins or CDS default:proteins

NR注释

# 优先级1 NCBI 和ensembl 等公共数据库上 下载 
# 优先级2 比对本物种 
# 优先级3 比对大库
Usage:
        /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/NR.sh -i cds.fa -o outdir -t CDS -s animal
OPTIONS:
        -h: help information
        -i: fasta file
        -o: output dir
        -t: proteins or CDS default:proteins
        -s: species [ animal plant fungi prokaryotes virus taxid(9606 10096 ..) ]

PPI注释

# 优先级1：脚本内会判断是否匹配上已知数据库
Usage:
        /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/NR.sh -i cds.fa -o outdir -t CDS -s animal
OPTIONS:
        -h: help information
        -i: fasta file
        -o: output dir
        -t: proteins or CDS default:proteins
        -s: species [ animal plant fungi prokaryotes virus taxid(9606 10096 ..) ]

转录因子注释

# 优先级1：脚本内会判断是否匹配上已知数据库
Usage:
        /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/NR.sh -i cds.fa -o outdir -t CDS -s animal
OPTIONS:
        -h: help information
        -i: fasta file
        -o: output dir
        -t: proteins or CDS default:proteins
        -s: species [ animal plant fungi prokaryotes virus taxid(9606 10096 ..) ]

合并注释文件

cd Annotation/tmp
perl /Business/psn_company/Work/Transcriptome/Datum/Public/Database/genome_prepare/script/merge_anno.v1.pl -g ../../GTF/Mus_musculus.GRCm39.111.gtf GO KEGG eggNOG eggNOG_Category Swissprot NR SwissprotName Entrez_geneID  BP CC MF Pathway -o ../Annotation.xls

# 统计：
perl /Business/psn_company/t01/public/ScriptManager/lib/PrepRef/stat_Annotation/stat_Annotation.pl -i Annotation.xls -o Annotation.stat

ncbi 参考基因组下载

说明：下载使用的是Aspera软件；设置的速度上限在100M；下载速度大概在 90M左右；下载完成后会整理好整理基因组流程需要的配置文件

如果有文件下载失败；则重新运行下载脚本；会跳过已下载的文件

/personalbio1/R_Report/Users/hudb/script/download_genome/ncbi.sh 
Usage:
	/personalbio1/R_Report/Users/hudb/script/download_genome/ncbi.sh -a [GCF.. GCA..] -o outdir 
assembly information: 245: /personalbio1/R_Report/Users/hudb/script/download_genome/ncbi/assembly_summary*txt
OPTIONS:
	-h: help information
	-a: assembly_accession
	-o: output dir

参数说明

-a NCBI assembly_accession编号；对应唯一的基因组；详细信息见 /personalbio1/R_Report/Users/hudb/script/download_genome/ncbi/assembly_summary*txt
-o 输出目录；

下载逻辑

。。。

ensembl 参考基因组下载

说明：ensembl 数据均使用wget下载；下载均加了-c参数；如果日志里有某些文件下载失败；则重新运行整个下载脚本

Usage: /personalbio1/R_Report/Users/hudb/script/download_genome/ensembl.sh -n homo_sapiens -v 111 -o ./ -r 3

Description:
Download ensembl genome;for all species to see /personalbio1/R_Report/Users/hudb/script/download_genome/ensembl/all_species.txt
Options:
-h, --help            Show this help message and exit.
-v, --version         ensembl genome version;default: current version
-o, --outdir          Output directory;default: ./
-n, --name            species name 
-r, --retry           curl and wget retry times;default: 3

参数说明

-v ensembl数据库基因组版本;默认数据库当前版本
-o 输出目录；
-n 物种名；是数据库ftp链接下的物种名；详细见/personalbio1/R_Report/Users/hudb/script/download_genome/ensembl/all_species.txt
-r 下载和请求网页的次数

Readme.md

真核有参参考基因组整理

非常规gff纠错和矫正为ensembl格式gtf

gff 转为 gtf

方法一：使用工具AGAT

方法二：手动处理；具体问题具体分析

参考基因组整理自动化流程

参数说明

配置文件 genome_conf.ini 介绍

参考基因组整理分步流程

参考基因组需要整理的文件 refCfg.ini 说明

DNA相关文件

GTF相关文件

annoDB文件夹

基因功能注释

GO注释

KEGG注释

eggnog注释

Swissprot注释

NR注释

PPI注释

转录因子注释

合并注释文件

ncbi 参考基因组下载

参数说明

下载逻辑

ensembl 参考基因组下载

参数说明

下载逻辑

只有基因组序列文件和基因结构注释gff文件

其他数据库下载

Readme.md

真核有参参考基因组整理

非常规gff纠错和矫正为ensembl格式gtf

gff 转为 gtf

方法一：使用工具AGAT

方法二：手动处理；具体问题具体分析

参考基因组整理自动化流程

参数说明

配置文件 genome_conf.ini 介绍

参考基因组整理分步流程

参考基因组需要整理的文件 refCfg.ini 说明

DNA相关文件

GTF相关文件

annoDB文件夹

基因功能注释

GO注释

KEGG注释

eggnog注释

Swissprot注释

NR注释

PPI注释

转录因子注释

合并注释文件

ncbi 参考基因组下载

参数说明

下载逻辑

ensembl 参考基因组下载

参数说明

下载逻辑

只有基因组序列文件和基因结构注释gff文件

其他数据库 下载

其他数据库下载