STAR输出的Aligned.toTranscriptome.out.bam深度计算

[trista1115]

您好， 我在使用STAR输出的Aligned.toTranscriptome.out.bam计算覆盖度和深度的时候，会出现下面的warning。想问pandepth在这种情况下是否可以使用？STAR用的也是这个gtf文件。

pandepth -t 3 -g gencode.vM32.basic.annotation.gtf -i Aligned.toTranscriptome.out.bam -o file

log文件显示： Warning: This region may be incorrect Warning: Can't find index file of input BAM/CRAM. PanDepth will run in No Index mode: Aligned.toTranscriptome.out.bam INFO: Input data read done

非常期待您的回复！谢谢！ 祝好

[HuiyangYu]

您好，非常抱歉给您造成的困扰。不知您是否可以提供gencode.vM32.basic.annotation.gtf的下载地址，我好进行检查。

[trista1115]

感谢您的快速回复，这是连接https://www.gencodegenes.org/mouse/release_M32.html gencode.vM32.basic.annotation.gtf对应的是这个gtf

[HuiyangYu]

出现这个warning的原因有2种： （1）gtf中的某个或者某些染色体在bam的header（参考基因组中的染色体名称）中无法找到； （2）gtf中某些区段的起始位置大于终止位置，我已经在gencode.vM32.basic.annotation.gtf进行了检查，这种情况是不存在的。 可能需要麻烦您通过samtools view -H Aligned.toTranscriptome.out.bam来检查bam中的染色体名称和数量和gtf文件中是否完全一致。

[HuiyangYu]

我猜测，核心的原因可能是比对的bam文件的染色体名称是没有'chr'的前缀的，也就是参考基因组的染色体名称为1，2，3 ...... M，X，Y，不知道是不是这个原因？

[trista1115]

感谢您的快速回复，我检查了Aligned.toTranscriptome.out.bam 的名称不是染色体，而是转录本id。可能需要我调整一下gtf文件的格式。 @SQ SN:ENSMUST00000193812.2 LN:1070 @SQ SN:ENSMUST00000082908.3 LN:110 @SQ SN:ENSMUST00000070533.5 LN:3634

[trista1115]

另外，虽然报了warning，但其实pandepth还是会计算出深度和覆盖度，那么这些结果是否属于不可信的结果？

[HuiyangYu]

您可以把使用STAR比对产生bam的命令，使用的参考基因组链接，pandepth计算的结果文件（随机重命名）贴（传）上来，我进行全面的检查。

[trista1115]

再次感谢您的回复！

star=2.7.3a

参考基因组的地址：https://www.gencodegenes.org/mouse/release_M32.html

fa_path="GRCm39.primary_assembly.genome.fa" gtf_path="gencode.vM32.basic.annotation.gtf" star_ref="STAR_genecode_m32" mkdir -p $star_ref

index

AA.chr.stat.gz STAR --runMode genomeGenerate \ --runThreadN 20 \ --genomeDir ${star_ref} \ --genomeFastaFiles ${fa_path} \ --sjdbGTFfile ${gtf_path} \ --sjdbOverhang 100

align

STAR --genomeDir ${star_ref} --readFilesIn AA_R1_val_1.fq.gz AA_R2_val_2.fq.gz --outFileNamePrefix AA \ --runThreadN 10 --readFilesCommand zcat --outSAMunmapped Within --outFilterType BySJout \ --outSAMattributes NH HI AS NM MD --outFilterMultimapNmax 20 --outFilterMismatchNmax 999 \ --outFilterMismatchNoverReadLmax 0.04 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 \ --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --sjdbScore 1 --genomeLoad LoadAndRemove --limitBAMsortRAM 10000000000 \ --outSAMtype BAM SortedByCoordinate --quantMode TranscriptomeSAM

pandepth

zcat AA.chr.stat.gz | head -n 5

Chr Length CoveredSite TotalDepth Coverage(%) MeanDepth

ENSMUST00000193812.2 1070 0 0 0.00 0.00 ENSMUST00000082908.3 110 0 0 0.00 0.00 ENSMUST00000070533.5 3634 836 1406 23.00 0.39 ENSMUST00000192857.2 480 147 294 30.62 0.61

[HuiyangYu]

非常抱歉，由于我自己暂时没有找到小鼠的数据，并且之前STAR比对的bam也找不到了，可能还需要麻烦您截下一部分bam（50-100M左右），并且上传上来。太不好意思了。

[HuiyangYu]

也可以直接截取最短的常染色体，截取的bam需要包含完整的header信息。

[trista1115]

不好意思，回复晚了，您看一下我上传的bam文件是否可用？我只截取了前面500个比对信息。 AA.500.bam.zip

[HuiyangYu]

我重新检查了你传的AA.500.bam/AA.chr.stat.gz，在bam文件中，所有的reads都是比对到了转录本ID中，此时你通过pandepth -i AA.bam -o AA得到的AA.chr.stat.gz就是对bam中的第3列进行统计。针对该bam文件，结果是完全正确的。但是由于不确定你为什么比对到转录本，而不是基因组，所以不确定这个统计结果对你后续的需求是否是合理的。

[trista1115]

正好借此机会跟您讨论一下， 主要是考虑到用的count table 是由rsem基于转录本的比对信息产生的，因此觉得reads比对到转录本上的覆盖度和深度会更好反应出数据情况？

[HuiyangYu]

这个问题可能几句话探讨不清楚，可以通过邮件或者其他方式来讨论。在这里可能不是非常合适。我邮件地址为：yuhuiyang1234@126.com。

[trista1115]

非常感谢您的快速回复！ 祝好！

[HuiyangYu]

把你的问题描述清楚，给我直接发邮件即可。

[hewm2008]

会洋 可以在readme里面把我的QQ 群放上 ，我设多一个管理员