SYSU-zhanglab / RNA-m5C

m5C mapping and site calling pipeline.
MIT License
8 stars 7 forks source link

关于fastq数据的疑问 #2

Closed fork16 closed 1 year ago

fork16 commented 1 year ago

刘博您好!非常感谢您提供分析的代码,让我非常清晰的理解了其中的原理,但我有一个问题还想请您指导一下:

我从GSE 下载的 m5C 数据内我发现一个规律,_1.fastq.gz 文件内的 G 碱基特别少,而 A 碱基特别多;_2.fastq.gz 文件内的 C 碱基特别少,而 T 碱基特别多。并且其他RNA m5C 的数据也是如此,所以我可以认为 1 端文件是 GA 转换, 2 端文件是 CT 转换吗?那我在运行hisat2 比对时 -F fwd.fastq 是否应该要输入 2端文件?因为二端是 CT 转换,正好对应Python 脚本内的 C2T。亦或是我的理解其实是错的,-F 应该输入1端文件?

恳请刘博能解答我的疑惑,非常感谢。

fork16 commented 1 year ago

还有一个问题是在 pdf 里面 cutadapt 有这句话 #Suppose using gunzip files as the input, dUTP stranded library (read1-ATC; read2-ATG) and standard illumina adapters ,这个括号内的 read1-ATC; read2-ATG 指的是什么呢?

jhfoxliu commented 1 year ago

这个与建库的方法有关系。dUTP-based 的directional RNA-Seq的read1是互补链,所以C多G少(C的互补)。如果用其他方法比如small RNA-seq,read1会与RNA链一致(G多C少)。

-F 指的是与RNA链一致的read。(如果只测了单端其read1是反向互补,需要先对read1的序列进行反向互补,质量值倒置)。

fork16 commented 1 year ago

非常感谢解答,已经了解了,所以ATC 指的是C多,ATG指的是G多的意思。

这里顺便提一下我遇到的一些问题及解决办法,供后来人参考:

在#Lift-over to genome (bowtie2_trans.bam to bowtie2_genome.bam)步骤中,这里会出现两个问题:

第一个是使用 Homo_sapiens.GRCh38.91.gtf(我用的是hg38的91版本)生成的 Homo_sapiens.GRCh38.91.anno 文件内的转录本ID是没有版本号的,而 Homo_sapiens.GRCh38.91.all.RNA.fa 文件内的转录本ID 是有版本号的,所以在这一步的时候会导致ID无法匹配从而导致 lift over 失败。解决办法是 把 fasta文件内ID 的版本号删除即可。

第二个是运行时会输出这个问题 [E::idx_find_and_load] Could not retrieve index file for 'bowtie2_trans.sorted.2.bam' ,虽然不会失败,但还是有解决办法。解决办法就是 使用 samtools sort 先进行排序,然后再用 samtools index 构建索引即可。

最后再次感谢刘博的帮助。

ponyaaaa commented 9 months ago

第二个问题我一开始也有碰到,之后改了pysam版本改到和作者pipline里一样之后就没有了;很感谢你的分享,刚好碰到lift over 失败这个问题!