三代转录组几大法宝:全长转录本、鉴定异构体(Isoforms)、鉴定可变剪接、鉴定基因融合,鉴定变异(SNP和INDEL),长读长带来二代测序不可比拟的优势!
第三代测序的常用的两大平台(Pacific Biosciences)、(Oxford Nanopore)。更为常用的是pacbio平台的实时单分子测序技术。对于融合基因的鉴定比二代更为准确。在植物领域中,融合基因的研究较少,融合基因定量软件体系也不成熟。该文章记录我对小麦的融合基因的定量。使用二代测序数据对三代转录组进行定量。
首先要得到由cDNA_Cupcake软件中的fusion_finderpy( cDNA_Cupcake链接 )脚本对sample_clusteredhqfasta文件(由pacbio官方提供的三代质控pipeline得到 IsoSeq链接 ),得到下图中的文件。
其中ALL_SAMPLErepfa文件的内容如下:
ALL_SAMPLE_correctedfasta文件内容如下,后面我们将从此文件提取两个基因片段的位置,从而确定融合连接点。>PBfusion11表示融合基因的左边序列,>PBfusion12表示融合基因的右边序列,根据左边序列长度从而确定junction位置。
1、构建索引文件,以ALL_SAMPLErepfa文件中的序列作为reference建立索引
--bowtie2 :调用 bowtie2软件作为比对软件,以构建索引文件
ALL_SAMPLErepfa:构建索引对象文件
<path>/name:<path>是储存索引文件路径,name是索引名称
结果如下:
2、RSEM开始比对定量
-p:线程数
--bowtie2:调用bowtie2软件进行比对
--bowtie2-mismatch-rate:比对错配率,默认为01
--paired-end:输入文件为二代双端测序文件,sample_R1fastpfqgz为二代测序上游文件,sample_R2fastpfqgz为二代测序下游文件
<path>/Wheat_fusion_ref:bowtie2构建索引文件位置及名称
<path_aim>/sample:生成结果的路径与文件前缀名
结果如下:
samplegenesresults:基因表达量
sampleisoformsresults:异构体表达量
sampletranscriptbam:比对结果bam文件
我们只需要得到有效比对的reads数量,而上面的比对结果文件以及定量文件都包含了无效比对和有效比对两种情况。下面对上面的结果文件进行提取有效比对结果。
1、首先将bam文件转为sam文件
2、筛选比对结果,使用下面python脚本
使用下面脚本需要:ALL_SAMPLE_correctedfasta、sampleisoformsresults文件
usage:python Exac_fusion_junction_sqpy /sample
上面脚本运行完成后会得到两个文件。
1、sampletranscrip_30alignedtxt 文件
2、sample30Isoformsresults
得到这两个文件,融合基因的表达定量就完成了。die "usage: $0 <fasta> [<fasta2> ] 1>seqfa 2>seqlen\n" if (@ARGV < 1);
die 表示终止脚本运行,并显示出die后面的双引号里面的内容。
die "" if (@ARGV < 1);则表示如果脚本运行时后面跟的参数少于1个,那么就停止运行并输出信息。
从你的这句来看,意思应该是你的脚本假设为runpl
那么运行的时候 在命令行输入的格式应该是 perl runpl <第一个FA文件名> <第二个FA文件名> 1>seqfa 2>seqlen
其中脚本后面跟的FA文件个数应该不限制,只要跟就可以了,可以1个也可以10个20个不限制。
然后将序列和长度分别存放在 seqfa 和seqlen两个文件里面。
其中1> 表示的是普通回显信息存放。
2>一般是用来将屏显错误的信息重定向到某个文件。 这个应该是这个脚本利用了这点输出的。
所以你也不用过多纠结了,只要明白 1>seqfa 表示把序列输出到seqfa 文件 2>seqlen表示把序列长度输出到seqlen就可以了。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)