生信练习题一

   一.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覆盖度

  二.

1）编写脚本或选择适当工具，统计 vcf 中变异位点的 Qual 值分布情况，并画图展示

zcat chr17.vcf.gz |grep -v '^#' | awk -F "	" '{if (FNR == 1) print "POS	qualq" ;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")

2）提取TP53基因的变异，说明变异位点是的数目（纯和和杂合）。

cat tp53_stat.sh
#! usr/bin/bash
for i in {10..12};do
#echo $i
cat tp54.vcf | grep '#CHROM' |awk -F "	" '{print $'$i'}'
echo -e "hom	 `cat tp54.vcf | grep '^chr17' |awk -F "	" '{print $'$i'}' | grep -v '0/0' | awk -F ":" '{print $1}' | awk -F "/" '{if($1==$2){print $0}}' | wc -l`"
echo -e "het	 `cat tp54.vcf | grep '^chr17' |awk -F "	" '{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 = '	',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	estfq_dataabc.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('
').split(' ')[0][1:]
                seq[seq_name]=''
            else:
                seq[seq_name]+=line.strip('
')
    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' + '	' + 'length' + '
')
    for k, v in seqname_length_sort.items():
        f.write(k + '	' + str(v) + '
')
    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 + "
" + "n50序列长度 : %s" % (seqname_length_sort[seqname]))
            break

calculate_n50()

结果：

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' + '	' + 'GC_percent' + '
')
    for k, v in location.items():
        f.write(k + '	' + str(v) + '
')
    print('over')
    f.close()

R画图展示

setwd("D:/pycharm_ngs_programs/hg19_stat")

d2 <- read.table(file = 'seqname_length_da2xiao.txt',sep = '	',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 = '	',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')

结果：

或者：

d1 <- read.table(file = 'Interval_GC_percent.txt',sep = '	',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)