Molecular Cell 文章 ribosome pausing 结果复现 (二) (PCR 去重)

人是一团欲望,满足不了便会难受

1引言

上一篇在去除通用接头的代码中,设置的 -m 参数指定最短保留的 reads 长度为 25nt, 但是文章里写的不是 15nt 吗?这里就需要注意我们去除通用接头的时候,实际上 reads 长度包括了(6+4)nt 的随机序列,所以 15+10 即为保留的最短长度。

前面我们去除了 umi 序列并加到了每条 read 的名称里,比对到基因组上后得到了 sam 文件, 今天对 sam 文件进行 PCR 去重。

2PCR 去重

先看看 sam 文件:

@SQ     SN:JH584295.1   LN:1976
@PG     ID:hisat2       PN:hisat2       VN:2.2.1        CL:"/home/zhoulab/anaconda3/bin/hisat2-align-s --wrapper basic-0 -p 20 -x 0.index-data/grcm38/genome>
SRR12594201.4|TAGGTNCGA 4       *       0       0       *       *       0       0       GACGAAGCCGAGCGCACGGGGTCGGCGGCGAT        EEEEEEEEEEEEEEEEEEEEEEEEEEEE>
SRR12594201.1|TAGAGNTNC 4       *       0       0       *       *       0       0       TCGACGAAGCCGAGCGCACGNGNTCGGCGGCGAT      /AEEEEEA6EA<EEAAE/E6#/#EEAEE>
SRR12594201.3|TAGCGNATC 4       *       0       0       *       *       0       0       TCGACGAAGCCGAGCGCACGGGGTCGGCGGC EEEEEAEE6EA<<AE/E//EE6EEE/EE/EE YT:Z

定义函数:julia语言编写

using XAM

# define function
function removePCRDup(inputFile,outputFile;outputType="txt")
    readInfoDict = Dict{String,String}()

    # load sam file
    reader = open(SAM.Reader, inputFile)

    # get sam header
    head = SAM.header(reader)
    # output to sam file
    samw = SAM.Writer(open(outputFile, "w"), head)

    record = SAM.Record()
    while !eof(reader)
        empty!(record)
        read!(reader, record)
        # do something
        if SAM.ismapped(record)
            # tags
            samName,randmer = split(SAM.tempname(record),"|")
            flag = SAM.flag(record)
            refNmae = SAM.refname(record)
            alignPos = SAM.position(record)
            readLength = SAM.seqlength(record)

            # not include duplicates
            key = join([randmer,refNmae,alignPos],"|")
            if !haskey(readInfoDict,key)
                readInfoDict[key] = ""

                # output data
                if outputType == "txt"
                    write(output,join([flag,refNmae,alignPos,readLength],"\t")*"\n")
                elseif outputType == "sam"
                    write(samw,record)
                else
                    println("Please give coorect type!")
                end
            else
                continue
            end
        end
    end

    # close file
    close(samw)
end

参数说明:

这里给了 outputType 可以是 txt 或者 sam 文件,如果输出 sam 文件可能文件会大一些,txt 格式只保留了几列重要的信息,文件会比较小,但是对于后面的分析足够了。

批量去重:

# sample
sample = ["SRR12594201","SRR12594202","SRR12594203","SRR12594204","SRR12594205","SRR12594206",
        "SRR12594207","SRR12594208","SRR12594209","SRR12594210","SRR12594211","SRR12594212",
        "SRR12594213","SRR12594214","SRR12594215","SRR12594216","SRR12594217","SRR12594218"]

for i in 1:length(sample)
    removePCRDup(join(["3.map-data/",sample[i],".sam"],""),
                 join(["3.map-data/",sample[i],".uniq.txt"],""),"txt")
end

txt 结果:

几列为 flag值,染色体,比对的位置,read长度:

16 4 114969279 29
16 16 61056439 29
0 5 146261023 31

我看了一个样本,大概去重后 reads 数量会减少一半左右,可以看到由于 pcr 扩增产生的重复还是挺多的。

3结尾

如果想输出 sam,输出格式修改一下即可,可以拿去重后的 sam 转 bam 去做其它分析或者可视化等。

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,数据代码已上传至QQ群,欢迎加入下载。

群二维码：

老俊俊微信：

知识星球：

所以今天你学习了吗？

今天的分享就到这里了，敬请期待下一篇！

最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！

如果觉得对您帮助很大，赏杯快乐水喝喝吧！

 往期回顾 

◀...