生信練習題一

時間 2019-12-01

原文原文鏈接

一. data/newBGIseq500_1.fq和data/newBGIseq500_2.fq中是基於BGIseq500測序平臺的一種真核生物基因組DNA的PE101測序數據，插入片斷長度爲450 bp；已知該基因組大小約在6M左右。

1）請統計本次測序的PE reads數是多少對reads？

2) 理論上可否使基因組99%以上的區域達到至少40X覆蓋？請簡要寫出推理和計算的過程與結果，數值計算使用R等工具時請寫出所用代碼。

wc -l data/newBGIseq500_1.fq|awk '{print $1/4}' 
#結果爲1599999

# 根據上一步結果計算所有的數據量 base<-1599999*200

#計算測序深度 dep<-base/6000000 #因爲基因組dna長度長，某片斷被檢測到的機率p<<1,而且測序過程當中會產生趨向無窮的reads，所以鹼基被測到的深度符合泊松分佈

#r語言的ppois()函數表示累積泊松分佈函數，所以，要計算基因組鹼基至少40X覆蓋的機率，應先求出0-39X被覆蓋的累積分佈機率ppois(39,dep),再用1-此機率就是大於等於40X的機率，

即： 1-ppois(39, dep) [1] 0.975145

#所以，理論上認爲，該數據量不能使基因組99%以上區域達到至少40X覆蓋度

二.python

1）編寫腳本或選擇適當工具，統計 vcf 中變異位點的 Qual 值分佈狀況，並畫圖展現bash

zcat chr17.vcf.gz |grep -v '^#' | awk -F "\t" '{if (FNR == 1) print "POS\tqualq" ;else print  $2,$6}' > chr17.qual

#! /usr/bin/env Rscript
##  Qual.R
png("qualstatistic-1.png")
qualddata <- read.table(file ="chr17.qual",sep = '',header = T) #注意時空格分隔
head(qualddata);dim(qualddata);class(qualddata)
colnames(qualddata) <- c('POS','QUAL')
p <- hist(qualddata$QUAL,main = "Qual Hist")
dev.off()

#方法二 ggplot
library(ggplot2)
ggplot(qualddata,aes(x=QUAL))+
# 直方圖函數：binwidth設置組距
geom_histogram(binwidth = 200000,fill='lightblue',color='black')+
xlab("QUAL")+ylab("Frequency")app

2）提取TP53基因的變異，說明變異位點是的數目（純和和雜合）。ide

cat tp53_stat.sh
#! usr/bin/bash
for i in {10..12};do
#echo $i
cat tp54.vcf | grep '#CHROM' |awk -F "\t" '{print $'$i'}'
echo -e "hom\t `cat tp54.vcf | grep '^chr17' |awk -F "\t" '{print $'$i'}' | grep -v '0/0' | awk -F ":" '{print $1}' | awk -F "/" '{if($1==$2){print $0}}' | wc -l`"
echo -e "het\t `cat tp54.vcf | grep '^chr17' |awk -F "\t" '{print $'$i'}' | grep -v '0/0' | awk -F ":" '{print $1}' | awk -F "/" '{if($1!=$2){print $0}}' | wc -l`"
done

bash tp53_stat.sh 
結果以下
27DMBDM4YT
hom     9
het     8
7XKZJA3JWX
hom     6
het     33
APRDKV0LDS
hom     8
het     8

cat tp53_select_out.txt
type    27DMBDM4YT    7XKZJA3JWX    APRDKV0LDS
hom    9    6    8
het    8    33    8

R畫圖
library(reshape)
library(ggplot2)

d1 <- read.table(file ="tp53_select_out.txt",sep = '\t',header = T)
head(d1)
#寬表格 整理成長表格
d1_m <-melt(d1,id.vars=c("type"))
head(d1_m)
ggplot(d1_m, aes(x=type, y=value)) +
geom_bar(stat="identity", position="dodge", aes(fill=variable))#position: dodge: 柱子並排放置


題三：
#下載並安裝 SOAPdenovo 軟件 ，設置 -K參數爲35對該數據 進行de novo組裝 ， #並畫出組裝結果序列從長到短的度累積曲線圖。 #計算組裝結果的 N50 #統計組裝結果GC 含量分佈，並畫圖展現

SOAPdenovo下載地址：
https://sourceforge.net/projects/soapdenovo2/files/latest/download
安裝：
make
軟件配置文件
#參考簡書https://www.jianshu.com/p/30c74297513a  SOAPdenovo 下載安裝 與使用  ，config 配置說明
config 配置說明:
除了輸入文件路徑外，還包含如下幾個參數的設置
1.    avg_ins
文庫插入片斷的平均長度，在實際設置時，能夠參考文庫size分佈圖，取峯值便可
2.     reverse_seq
是否須要將序列反向互補，對於pair-end數據，不須要反向互補，設置爲0；對於mate-pair數據，須要反向互補，設置爲1
3.    asm_flags
1表示只組裝contig. 2表示只組裝scaffold,3表示同時組裝contig和scaffold,4表示只補gap
4.    rd_len_cutof
序列長度閾值，做用和max_rd_len相同，大於該長度的序列會被切除到該長度
5.    rank
設置不一樣文庫數據的優先級順序，取值範圍爲整數，rank值相同的多個文庫，在組裝scaffold時，會同時使用。
6.    pair_num_cutoff
contig或者scaffold以前的最小overlap個數，對於pair-end數據，默認值爲3；對於mate-paird數據，默認值爲5
7.    map_len
比對長度的最小閾值，對於pair-end數據，默認值爲32；對於mate-pair數據，默認值爲35

配置示例以下：cat config_file

#cat config_file
max_rd_len=100
[LIB]
avg_ins=450
asm_flag=3
map_len=32
#####    rd_len_cutoff=100 序列長度閾值，做用和max_rd_len相同，大於該長度的序列會被切除到該長度
pair_num_cufoff=3
reverse_seq=0
rank=1
q1=/home_extend/u****/R/exam/newBGIseq500_1.fq   #注意使用上一步過濾後的 clean.fq
q2=/home_extend/u****/R/exam/newBGIseq500_2.fq

運行
./SOAPdenovo-63mer all -s config_file -K 35 -R -o /home_extend/u1010/R/exam/abc

結果函數

Version 2.04: released on July 13th, 2012
Compile Jul  9 2013     11:57:16

********************
Pregraph
********************

Parameters: pregraph -s config_file -K 35 -R -o /home_extend/u****/R/exam/abc

In config_file, 1 lib(s), maximum read length 100, maximum name length 256.

8 thread(s) initialized.
.......

累計曲線與N50 計算：工具

#!/usr/bin/env/python
#encoding=utf8
import sys
#args = sys.argv
#fasta_file = args[1]
fasta_file = r'E:\pylab\testfq_data\abc.scafSeq'
def read_fasta():
    with open(fasta_file) as F_in:
        seq_name= ''
        seq = {}
        for line in F_in:
            if line.startswith(">"):
                seq_name=line.strip('\n').split(' ')[0][1:]
                seq[seq_name]=''
            else:
                seq[seq_name]+=line.strip('\n')
    return seq

def get_seqname_len():
    seq = read_fasta()
    seqname_length = {}
    seqname_length_sort ={}
    for k, v in seq.items():
        seqname_length[k] = len(v)
    seqname_length_sort = sorted(seqname_length.items(), key=lambda x: x[1])#排序
    seqname_length_sort = dict(seqname_length_sort)
    return seqname_length_sort

def write_seqname_length():
    seqname_length_sort = get_seqname_len()
    f = open('t5.txt', 'w')
    f.write('seq_name' + '\t' + 'length' + '\n')
    for k, v in seqname_length_sort.items():
        f.write(k + '\t' + str(v) + '\n')
    f.close()


def calculate_n50():
    sum_len = 0
    n50 = 0
    seqname_length_sort = get_seqname_len()
    for seqname,length in seqname_length_sort.items():
        sum_len +=length
        L.append(length)
    for seqname, length in seqname_length_sort.items():
        n50 += length
        count +=1
        if n50 >= sum_len/2:
            print("n50序列名 : %s" % seqname + "\n" + "n50序列長度 : %s" % (seqname_length_sort[seqname]))
            break

calculate_n50()

結果：spa

n50序列名 : scaffold204
n50序列長度 : 40176

根據t5.txt 文件內容可知
#       seq_name        length
# min =  C156737           100
# max = scaffold82      176232

統計該區間的GC percent,存到文件 Interval_GC_percent.txt 中R做圖

seq = read_fasta()
def Interval_GC_count():
    min = 100
    max = 176232
    extend = 1000
    location = {}
    while min < max :
    　　 count = 0
        SUM_GC = 0
        SUM_LEN = 0
        regin = range(min, min + 1000)
        for name ,sequence in seq.items():
            length = len(sequence)
            GC  = sequence.count('G') + sequence.count('C')
            if length in regin:
                #count +=1
                SUM_GC += GC
                SUM_LEN += length
        if SUM_LEN >0:
            OUT ="%.2f"%(SUM_GC/SUM_LEN)
            location[str(min) + '-' + str(min + 1000)] = OUT
        min +=1000
        if min > max:
            break
    return location
location= Interval_GC_count()
def WRITE():
    f = open('Interval_GC_percent.txt', 'w')
    f.write('Interval' + '\t' + 'GC_percent' + '\n')
    for k, v in location.items():
        f.write(k + '\t' + str(v) + '\n')
    print('over')
    f.close()

R畫圖展現.net

setwd("D:/pycharm_ngs_programs/hg19_stat")

d2 <- read.table(file = 'seqname_length_da2xiao.txt',sep = '\t',header = F)
head(d2)
cumsum_data <- cumsum(d2[,2]) #cumsum累加數據
cumsumframe <- as.data.frame(cumsum_data) #轉換成數據框
head(cumsumframe);dim(cumsumframe);class(cumsumframe)
df <- cbind.data.frame(cumsumframe,d2)#合併數據框

head(df);dim(df);class(df)
plot(df$cumsum_data,type="l",ylab='Total length',xlab = 'Seq num',main="序列長度累積曲線圖")

#或者
l<-seq(1,nrow(df),by=1 )
plot(x=l,y=df$cumsum_data,type="l",ylab='Total length',xlab = 'num',main="序列長度累積曲線圖")

readslength 區間與GC_percent

d1 <- read.table(file = 'Interval_GC_percent.txt',sep = '\t',header = T)
head(d1)
order_level<- d1$Interval
d1$Interval<- factor(d1$Interval,levels = order_level,ordered = T)
str(d1)
library(ggplot2)
ggplot(d1,aes(x=d1$Interval,y=GC_percent))+
  geom_bar(stat="identity",fill = "blue")+
  theme(axis.text.x=element_text(angle=90,hjust = 1,vjust = 1))+
  xlab("Interval(bp)")+ylab("GC_percent(%)")
  ggsave('Interval_gc_percent.png')

結果：code

或者：blog

d1 <- read.table(file = 'Interval_GC_percent.txt',sep = '\t',header = T)
head(d1);dim(d1)
inter <- d1$Interval
GC_percent <- d1$GC_percent
plot(GC_percent,type="b",ylim=c(0.2,0.6),xaxt="n",yaxt="n",bty ='o',
     main="bty默認取\"o\"",xlab="長度",ylab="GC_percent")
axis(side = 1,at=1:71,labels=inter,tick=TRUE,las = 2)