cell 代码优化测试复现一
测试开头




测试结尾

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


1前言
于光明中寻找黑暗,于黑暗中追寻光辉。
2引言
关于前面 cell 文章的复现,作者的数据是 比对到基因组上 的,此外 代码也有部分问题, 这里我采用 比对到转录组序列上 进行分析,对结果进行重现,然后优化修正作者的部分代码,再使其适用于转录组序列上的数据
。
往期复现:
3转录组序列准备
因为酵母多为
单个连续的 CDS 序列
,且大部分没有 UTR 序列
,为了方便分析,对其 CDS 序列两端分别拓展 50nt 长的序列作为转录组序列。
序列获取参考往期内容:
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)

7计算 read density
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!