search

MGI-tech-bioinformatics / DNBelab_C_Series_HT_scRNA-analysis-software

An open source and flexible pipeline to analysis high-throughput DNBelab C Series single-cell RNA datasets
MIT License
65 stars 24 forks source link

How to add a marker gene for building a custom reference for dnbc4tools? #85

Closed jiangjasson closed 3 months ago

jiangjasson commented 3 months ago

你好，我在构建人的带荧光参考基因组的时候遇到了个问题，我使用的是最新的Gencode的fasta和gtf文件，在run的时候发现添加的一个荧光基因会被跳过，rna run输出的信息如下：

2024-06-27 17:13:42,165 - data - INFO - Executing command: /home/dell/anaconda3/envs/dnbc4tools/lib/python3.8/site-packages/dnbc4tools/software/Anno -I /home/dell/Documents/jiangss/project/H1-Ascl1-mCherry-clone4/H1ASCL1PNEC/01.data/Aligned.out.bam -a /home/dell/Documents/jiangss/reference/GencodeRef/DNBC4_format/gencode.v46.primary_assembly.annotation.filter.mCherry.gtf -L /home/dell/Documents/jiangss/project/H1-Ascl1-mCherry-clone4/H1ASCL1PNEC/01.data/cDNA.barcode.counts.txt -o /home/dell/Documents/jiangss/project/H1-Ascl1-mCherry-clone4/H1ASCL1PNEC/01.data -c 90 -m chrM -B /home/dell/anaconda3/envs/dnbc4tools/lib/python3.8/site-packages/dnbc4tools/config/cellbarcode/scRNA_beads_darkReaction.json --anno 1 --intron
2024-06-27 18:54:37,411 - data - INFO - 2024-06-27 17:14:06.789 W process: Multiple gene IDs for gene mCherry: mCherry -- skipping

还有另外两个也被跳过了：

2024-06-27 17:14:06.851 W process: Strand disagreement for gene LINC02256 -- skipping
2024-06-27 17:14:08.156 W process: Strand disagreement for gene LINC03025 -- skipping

这是log目录的文件：20240627.txt

我添加了荧光基因的gtf文件如下：

(base) dell@dell-Precision-7920-Tower:~/Documents/jiangss/reference/GencodeRef/DNBC4_format$ tail -n 2 gencode.v46.primary_assembly.annotation.filter.mCherry.gtf 
KI270734.1  ENSEMBL UTR 138082  138482  .   -   .   gene_id "ENSG00000277196.4"; transcript_id "ENST00000621424.4"; gene_type "protein_coding"; gene_name "ENSG00000277196"; transcript_type "protein_coding"; transcript_name "ENST00000621424"; exon_number 15; exon_id "ENSE00003753010.1"; level 3; protein_id "ENSP00000481127.1"; transcript_support_level "1"; tag "basic"; tag "Ensembl_canonical";
chrF    custom  exon    1   708 .   +   .   gene_id "mCherry"; transcript_id "mCherry"; gene_name "mCherry"; gene_type "protein_coding";

添加了荧光基因的fasta文件如下：

(base) dell@dell-Precision-7920-Tower:~/Documents/jiangss/reference/GencodeRef/DNBC4_format$ tail -n 13 GRCh38.primary_assembly.genome.mCherry.fa 
>mCherry
ATGGTGAGCAAGGGCGAGGAGGATAACATGGCCATCATCAAGGAGTTCATGCGCTTCAAG
GTGCACATGGAGGGCTCCGTGAACGGCCACGAGTTCGAGATCGAGGGCGAGGGCGAGGGC
CGCCCCTACGAGGGCACCCAGACCGCCAAGCTGAAGGTGACCAAGGGTGGCCCCCTGCCC
TTCGCCTGGGACATCCTGTCCCCTCAGTTCATGTACGGCTCCAAGGCCTACGTGAAGCAC
CCCGCCGACATCCCCGACTACTTGAAGCTGTCCTTCCCCGAGGGCTTCAAGTGGGAGCGC
GTGATGAACTTCGAGGACGGCGGCGTGGTGACCGTGACCCAGGACTCCTCCCTGCAGGAC
GGCGAGTTCATCTACAAGGTGAAGCTGCGCGGCACCAACTTCCCCTCCGACGGCCCCGTA
ATGCAGAAGAAGACCATGGGCTGGGAGGCCTCCTCCGAGCGGATGTACCCCGAGGACGGC
GCCCTGAAGGGCGAGATCAAGCAGAGGCTGAAGCTGAAGGACGGCGGCCACTACGACGCT
GAGGTCAAGACCACCTACAAGGCCAAGAAGCCCGTGCAGCTGCCCGGCGCCTACAACGTC
AACATCAAGTTGGACATCACCTCCCACAACGAGGACTACACCATCGTGGAACAGTACGAA
CGCGCCGAGGGCCGCCACTCCACCGGCGGCATGGACGAGCTATATAAG

这个fasta荧光基因这里，我还试过>chrF FL，也是会被跳过。

请问是我哪里设置出问题了吗？

lishuangshuang0616 commented 3 months ago

To add entries to the GTF, you need at least gene/transcript, and exon information. The lack of gene information causes it to be recognized within the gene ENSG00000277196. Please adjust accordingly, ensuring to use tab separators. Also, the FASTA name should be chrF.

chrF    custom  gene    1   708 .   +   .   gene_id "mCherry"; gene_name "mCherry"; gene_type "protein_coding";
chrF    custom  exon    1   708 .   +   .   gene_id "mCherry"; transcript_id "mCherry"; gene_name "mCherry"; gene_type "protein_coding";
>chrF
ATGGTGAGCAAGGGCGAGGAGGATAACATGGCCATCATCAAGGAGTTCATGCGCTTCAAG
GTGCACATGGAGGGCTCCGTGAACGGCCACGAGTTCGAGATCGAGGGCGAGGGCGAGGGC
CGCCCCTACGAGGGCACCCAGACCGCCAAGCTGAAGGTGACCAAGGGTGGCCCCCTGCCC
TTCGCCTGGGACATCCTGTCCCCTCAGTTCATGTACGGCTCCAAGGCCTACGTGAAGCAC
CCCGCCGACATCCCCGACTACTTGAAGCTGTCCTTCCCCGAGGGCTTCAAGTGGGAGCGC
GTGATGAACTTCGAGGACGGCGGCGTGGTGACCGTGACCCAGGACTCCTCCCTGCAGGAC
GGCGAGTTCATCTACAAGGTGAAGCTGCGCGGCACCAACTTCCCCTCCGACGGCCCCGTA
ATGCAGAAGAAGACCATGGGCTGGGAGGCCTCCTCCGAGCGGATGTACCCCGAGGACGGC
GCCCTGAAGGGCGAGATCAAGCAGAGGCTGAAGCTGAAGGACGGCGGCCACTACGACGCT
GAGGTCAAGACCACCTACAAGGCCAAGAAGCCCGTGCAGCTGCCCGGCGCCTACAACGTC
AACATCAAGTTGGACATCACCTCCCACAACGAGGACTACACCATCGTGGAACAGTACGAA
CGCGCCGAGGGCCGCCACTCCACCGGCGGCATGGACGAGCTATATAAG
jiangjasson commented 3 months ago

非常感谢，按照你的步骤，我的荧光基因没有被跳过了，最后也成功生成了表达矩阵，但是我查看了下output/filter_matrix的features.tsv.gz(features.tsv.gz)发现有ENSG00000xxxxxx，这是什么原因呢？

(base) dell@dell-Precision-7920-Tower:~/Documents/jiangss/project/H1-Ascl1-mCherry-clone4/H1ASCL1PNEC/output/filter_matrix$ zless -S features.tsv.gz | head
ENSG00000238009
MIR1302-2HG
ENSG00000290826
ENSG00000290385
ENSG00000291215
ENSG00000293331
ENSG00000241860
ENSG00000241599
ENSG00000292994
ENSG00000235146

而且有12529个

(base) dell@dell-Precision-7920-Tower:~/Documents/jiangss/project/H1-Ascl1-mCherry-clone4/H1ASCL1PNEC/output/filter_matrix$ zgrep "^ENSG00" features.tsv.gz |wc
  12529   12529  200464

我对比了下公司给我的表达矩阵的features.tsv.gz(features.tsv.gz)，公司给的是正常的

(base) dell@dell-Precision-7920-Tower:/media/dell/My Passport/H1PNEC$ zless -S features.tsv.gz | head
MIR1302-2HG
AL627309.1
AL627309.5
AL627309.4
AL732372.1
AC114498.1
LINC01409
LINC01128
AL645608.6
AL390719.3

我看了下公司生成的html报告，公司的Species显示的是Human，而我用的是Homo_sapiens

我不知道是不是因为我用的是最新的fasta和gtf导致的，还是因为mkref或者rna run有参数没设置正确，还请麻烦老师帮我看下。

这是我mkref的命令：

dnbc4tools rna mkref --ingtf /home/dell/Documents/jiangss/reference/GencodeRef/DNBC4_format/gencode.v46.primary_assembly.annotation.filter.mCherry.gtf --fasta /home/dell/Documents/jiangss/reference/GencodeRef/DNBC4_format/GRCh38.primary_assembly.genome.mCherry.fa --threads 90 --species Homo_sapiens --genomeDir GRCh38-gencodev46-mCherry708 1>dnbcrnamkref.log 2>&1

这是rna run的命令：

dnbc4tools rna run --cDNAfastq1 /mnt/disk1/dataset/H1-ASCL1-Clone4-PNEC/rawdata/cDNA/H1Ascl1clone4_S1_L001_1.fq.gz --cDNAfastq2 /mnt/disk1/dataset/H1-ASCL1-Clone4-PNEC/rawdata/cDNA/H1Ascl1clone4_S1_L001_2.fq.gz --oligofastq1 /mnt/disk1/dataset/H1-ASCL1-Clone4-PNEC/rawdata/oligo/H1Ascl1clone4_S1_L001_1.fq.gz --oligofastq2 /mnt/disk1/dataset/H1-ASCL1-Clone4-PNEC/rawdata/oligo/H1Ascl1clone4_S1_L001_2.fq.gz --genomeDir /home/dell/Documents/jiangss/reference/dnbc4-genome/GRCh38-gencodev46-mCherry708 --name H1ASCL1PNEC --threads 90 1>dncbrnarun.log 2>&1
lishuangshuang0616 commented 3 months ago

应该只是gtf版本不一致的问题，可以看下v46版本相比于quick start中推荐使用的v32增加了哪些信息。 dnbc4tools使用的是gene_nametag信息作为feature。可以看到这些ENS大多数为LncRNA。

jiangjasson commented 3 months ago

感谢老师这么晚还回复，确实是gtf版本的问题，换成v32后就没有ENS了