cell 代码优化测试复现一

测试开头

没关注？伸出手指点这里—

1前言

于光明中寻找黑暗,于黑暗中追寻光辉。

2引言

关于前面 cell 文章的复现,作者的数据是 比对到基因组上 的,此外 代码也有部分问题, 这里我采用 比对到转录组序列上 进行分析,对结果进行重现,然后优化修正作者的部分代码,再使其适用于转录组序列上的数据。

往期复现:

跟着 Cell 学 Ribo-seq 分析一

跟着 Cell 学 Ribo-seq 分析二

跟着 Cell 学 Ribo-seq 分析三 (Metagene Plot)

跟着 Cell 学 Ribo-seq 分析四 (终)

3转录组序列准备

因为酵母多为单个连续的 CDS 序列,且大部分没有 UTR 序列,为了方便分析,对其 CDS 序列两端分别拓展 50nt 长的序列作为转录组序列。

序列获取参考往期内容:

提取酵母两端扩展 50nt 的 CDS 序列

4建立索引

对转录组序列建立索引:

# 7.map to trancriptome
mkdir trans-index
mkdir 6.map-trans-data

bowtie-build extened-CDS-seq.fa trans-index/extened-CDS-seq

5比对

使用去除 rRNA 和 tRNA 的测序数据进行比对:

for i in SRR5188601 SRR5188602 SRR5188603 SRR5188604 SRR5188613 SRR5188614 SRR5188615 SRR5188616
do
    bowtie -p 20 ./index/trans-index/extened-CDS-seq 
           -m 1 -v 2 --best --strata 
           <(zcat 4.rmrtRNA-data/${i}.rmtRNA.fq.gz) 
           6.map-trans-data/${i}.map
done

6分析 read 长度分布

定义函数:

# get read length from .map data
def readlength(inputFile,outputFile):
    # open file
    outfile = open(outputFile,'w')
    len_dict = {}
    with open(inputFile,'r') as input:
        for line in input:
            fileds = line.split()
            length = len(fileds[5])
            len_dict.setdefault(length,0)
            len_dict[length] += 1

    # sort dict
    len_dict_sorted = list(len_dict.items())
    len_dict_sorted.sort()

    # output
    for elem in len_dict_sorted:
        outfile.write('t'.join([str(elem[0]),str(elem[1])]) + 'n')
    outfile.close()

批量提取:

import os

# make folder
os.mkdir('./1.read-length-data')

# readlength('ssb1-inter-rep1.map','test.length.txt')

sample = ['ssb1-inter-rep1.map','ssb1-inter-rep2.map','ssb2-inter-rep1.map','ssb2-inter-rep2.map',
          'ssb1-trans-rep1.map','ssb1-trans-rep2.map','ssb2-trans-rep1.map','ssb2-trans-rep2.map']

# run
for i in range(0,8):
    readlength(sample[i],''.join(['1.read-length-data/',sample[i],'.len.txt']))

R 可视化

library(tidyverse)
library(ggplot2)

file <- list.files('1.read-length-data/',pattern = '.txt')

# bacth read
map_df(file,function(x){
  tmp = read.table(paste('1.read-length-data/',x,sep = ''))
  colnames(tmp) <- c('readlength','numbers')
  tmp$percent <- tmp$numbers/sum(tmp$numbers)
  tmp$sample <- sapply(strsplit(x,split = '\.'),'[',1)
  # filter length 25-35 nt
  tmp <- tmp %>% filter(readlength >= 25 & readlength <= 35)
  return(tmp)
}) -> len

# plot 6x15
ggplot(len,aes(x = factor(readlength),y = percent)) +
  geom_col() +
  theme_bw(base_size = 14) +
  xlab('') +
  facet_wrap(~sample,scales = 'free',ncol = 4)