有参转录组数据分析,样品和基因组比对效率只有零点几,这是为什么?

老师,您好,根据咱们转录组课程进行了实操,质控等正常,但是转录组数据比对到基因组这个步骤总是出现比对率只有零点几的情况,尝试了好几个转录组数据,从来没正常过。

水稻转录组测序数据显示如下,NCBI下载的单端数据:1 @SRR22419667.1 A00133:246:HNTYVDSXY:1:1101:30734:1454 length=300

      2 TCCCGTTGAGGTTACACTCCAAATTGAATGATGAGGCAACGAATAACAAAGCAGCACTCAGAAATCCATGATGGTCCAGGCATCTACAAATGATATATCGCGCAACTCCATTTGAGTACAAATGGGTCCAAGCCCTGGCAACAGGTGTATCCTTAAGAG

      3 +SRR22419667.1 A00133:246:HNTYVDSXY:1:1101:30734:1454 length=300

      4 FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

      5 @SRR22419667.2 A00133:246:HNTYVDSXY:1:1101:6822:1689 length=300

      6 GAGCTCAGTAACAGCAGTCTCTGAGTCAGCATGCCATGTGAGTAAATCAACAAAAGGACAAGTACCAGGATCTAGCAAATTGGGTCAACTAGCTACACGGAGAAATTTAGAGACTGAAGTTTTGCAAAGTCAGCAAGTTAGCAGCCCAGCCAGTTCCAC

      7 +SRR22419667.2 A00133:246:HNTYVDSXY:1:1101:6822:1689 length=300

      8 FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF.

#水稻基因组的文件从ensemble plants上下载,构建索引的代码为

#基因组信息下载到服务器/home/us024/Imranrnaseq/my_rnaseq/ref

wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-61/fasta/oryza_sativa/dna/Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz

wget https://ftp.ensemblgenomes.ebi.ac.uk/pub/plants/release-61/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.61.gff3.gz

gunzip -c Oryza_sativa.IRGSP-1.0.dna.toplevel.fa.gz > Oryza_sativa.IRGSP-1.0.dna.toplevel.fa

gunzip -c Oryza_sativa.IRGSP-1.0.61.gff3.gz > Oryza_sativa.IRGSP-1.0.61.gff3

#修改gff3文件

sed 's#gene:##'  Oryza_sativa.IRGSP-1.0.61.gff3|ee

sed 's#gene:##'  Oryza_sativa.IRGSP-1.0.61.gff3|sed 's#transcript:##' |ee

#然后把修改后的另存为另一个名字,运行如下

sed 's#gene:##'  Oryza_sativa.IRGSP-1.0.61.gff3|sed 's#transcript:##' > Oryza_sativa.IRGSP-1.0.61.gff31

#把源文件覆盖掉

mv Oryza_sativa.IRGSP-1.0.61.gff31  Oryza_sativa.IRGSP-1.0.61.gff3

#构建索引

sh $scriptdir/index.sh Oryza_sativa.IRGSP-1.0.dna.toplevel.fa Oryza_sativa.IRGSP-1.0.61.gff3

#构建索引后,gtf的部分显示,其它的生成的结果看着也没什么问题

 1 1       RAP2022-09-01   transcript      2983    10815   .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

      2 1       RAP2022-09-01   exon    2983    3268    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

      3 1       RAP2022-09-01   exon    3354    3616    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

      4 1       RAP2022-09-01   exon    4357    4455    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

      5 1       RAP2022-09-01   exon    5457    5560    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

      6 1       RAP2022-09-01   exon    7136    7944    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

      7 1       RAP2022-09-01   exon    8028    8150    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

      8 1       RAP2022-09-01   exon    8232    8320    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

      9 1       RAP2022-09-01   exon    8408    8608    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

     10 1       RAP2022-09-01   exon    9210    9615    .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

     11 1       RAP2022-09-01   exon    10102   10187   .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

     12 1       RAP2022-09-01   exon    10274   10430   .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

     13 1       RAP2022-09-01   exon    10504   10815   .       +       .       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

     14 1       RAP2022-09-01   CDS     3449    3616    .       +       0       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

     15 1       RAP2022-09-01   CDS     4357    4455    .       +       0       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

     16 1       RAP2022-09-01   CDS     5457    5560    .       +       0       transcript_id "Os01t0100100-01"; gene_id "Os01g0100100";

#构建索引后,gff3的部分显示

 69 #!genebuild-last-updated 2019-06

     70 1       IRGSP-1.0       chromosome      1       43270923        .       .       .       ID=chromosome:1;Alias=Chr1,AP014957.1,NC_029256.1

     71 ###

     72 1       RAP2022-09-01   gene    2983    10815   .       +       .       ID=Os01g0100100;biotype=protein_coding;gene_id=Os01g0100100;logic_name=rapdb_genes

     73 1       RAP2022-09-01   mRNA    2983    10815   .       +       .       ID=Os01t0100100-01;Parent=Os01g0100100;biotype=protein_coding;tag=Ensembl_canonical;tra

     74 1       RAP2022-09-01   exon    2983    3268    .       +       .       Parent=Os01t0100100-01;Name=Os01t0100100-01-E1;constitutive=1;ensembl_end_phase=-1;ense

     75 1       RAP2022-09-01   five_prime_UTR  2983    3268    .       +       .       Parent=Os01t0100100-01

     76 1       RAP2022-09-01   five_prime_UTR  3354    3448    .       +       .       Parent=Os01t0100100-01

     77 1       RAP2022-09-01   exon    3354    3616    .       +       .       Parent=Os01t0100100-01;Name=Os01t0100100-01-E2;constitutive=1;ensembl_end_phase=0;ensem

     78 1       RAP2022-09-01   CDS     3449    3616    .       +       0       ID=CDS:Os01t0100100-01;Parent=Os01t0100100-01;protein_id=Os01t0100100-01

     79 1       RAP2022-09-01   exon    4357    4455    .       +       .       Parent=Os01t0100100-01;Name=Os01t0100100-01-E3;constitutive=1;ensembl_end_phase=0;ensem

     80 1       RAP2022-09-01   CDS     4357    4455    .       +       0       ID=CDS:Os01t0100100-01;Parent=Os01t0100100-01;protein_id=Os01t0100100-01

     81 1       RAP2022-09-01   exon    5457    5560    .       +       .       Parent=Os01t0100100-01;Name=Os01t0100100-01-E4;constitutive=1;ensembl_end_phase=2;ensem

 

#比对到基因组,先比对一个

hisat2 -p 3 --rg-id=400ppm_rep1 --rg SM:400ppm_rep1 --rg LB:400ppm_rep1 --rg PL:ILLUMINA -x $REF_INDEX --dta -U $workdir/2.data_qc/400ppm_rep1.clean.fq.gz -S 400ppm_rep1.sam 2> 400ppm_rep1.summary

#比对结果

24819002 reads; of these:

  24819002 (100.00%) were unpaired; of these:

    24775361 (99.82%) aligned 0 times

    39911 (0.16%) aligned exactly 1 time

    3730 (0.02%) aligned >1 times

0.18% overall alignment rate

请先 登录 后评论
  • 0 关注
  • 0 收藏,50 浏览
  • naliupeony 提出于 3天前

相似问题