首先得到每條序列的長度豺憔，在這里使用seqkit軟件证舟。
seqkit軟件是一個強(qiáng)大的序列處理工具释树，安裝方法參見官方網(wǎng)站.
代碼如下：
seqkit fx2tab -j 30 -l -n -i -H file.fastq.gz > Length.txt
#參數(shù)：
# -j 是線程數(shù)
#-B, --base-content value 要輸出的堿基含量e.g. -B AT -B N
#-g, --gc print GC content
#-l, --length print sequence length
#-n, --name only print names
#-i, --only-id print ID instead of full head
head Length.txt # 查看Length.txt
結(jié)果如下所示：
Length.txt
提取長度分布：
seqkit fx2tab -j 30 -l -n -i -H file.fastq.gz |cut -f 2 > length.txt
將得到結(jié)果length.txt導(dǎo)入到R中繪制序列長度分布圖茄袖，
代碼如下：
library(tidyverse)
length <- read_tsv("Length.txt") %>% group_by(length) %>% summarise(Count = n())
sum <- sum(length$Count)
ggplot(length) +
geom_col(aes(length, Count), width = 0.8) +
geom_line(aes(length, Count), group = 1) +
geom_point(aes(length, Count)) +
scale_y_continuous(sec.axis = sec_axis(~.*100/sum, name = "% Relative Abundance")) +
xlab("R1.Length") +
theme_bw() + theme(panel.grid = element_blank(),
axis.title = element_text(size = 15))
ggsave("Length.png", height = 5, width = 8)
ggsave("Length.pdf", height = 5, width = 8)
結(jié)果如下：
圖1
上圖由于數(shù)據(jù)太多桥狡，去掉點(diǎn)和連線更清晰
library(tidyverse)
length <- read_tsv("Length.txt") %>% group_by(length) %>% summarise(Count = n())
sum <- sum(length$Count)
ggplot(length) +
geom_col(aes(length, Count), width = 0.8) +
# geom_line(aes(length, Count), group = 1) +
# geom_point(aes(length, Count)) +
scale_y_continuous(sec.axis = sec_axis(~.*100/sum, name = "% Relative Abundance")) +
xlab("Length") +
theme_bw() + theme(panel.grid = element_blank(),
axis.title = element_text(size = 15))
ggsave("Length.png", height = 5, width = 8)
ggsave("Length.pdf", height = 5, width = 8)
結(jié)果如下：
圖2
---------------------------------------------------------------------------------------------------------------------------------------------------I`m a line ! Thanks !-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
