搜索
查看: 2212|回复: 0

[Other] Biostar:课程23、24

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-9-2 20:39:17 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:István Albert
原文链接:https://mp.weixin.qq.com/s/HPwonZfGu_97QluA8Le6Ew

翻译:生物女博士
注:本系列课程翻译已获得原作者授权。

#为作者的注释
#*为译者的注释


Lecture 23 - 从埃博拉数据中Call SNPs

# 从多个样品中Call SNPS
# 从埃博拉项目中获取多个数据集。
# Ouch! 数据是以另一个序列作为参考来比对的。
# 那我们准备一个新的参考序列吧,没别的解决方法。
#* NCBI的这套软件最近运行都有问题,建议直接到NCBI下载好再上传服务器。


efetch -id KM034562 -format gb -db nucleotide > KM034562.gb
readseq -format=FASTA -o ~/refs/852/KM034562.fa
KM034562.gb
readseq -format=GFF -o KM034562.gff KM034562.gb
bwa index ~/refs/852/KM034562.fa
samtools faidx ~/refs/852/KM034562.fa



# 获取运行信息文件
esearch -db sra -query PRJNA257197  | efetch -format runinfo > runinfo.csv

# 选一些SRR(如果觉得冒险的话,可以弄多一些)
cat runinfo.csv | cut -f 1 -d ',' | grep SRR | head -10 > sraids.txt

# 我们把数据放到另一个文件夹中。
# 将会有很多文件。
mkdir -p sra
mkdir -p bam

# 从文件读取每个SRA ID
cat sraids.txt | awk ' { print "fastq-dump -X 20000 --split-files -O sra " $1 } ' > get-data.sh
#* 查看一下get-data.sh的内容
# 这个脚本将会执行数据
#* 数据的下载,以及将格式从sra变为fastq。
#* 解释一下这条命令的意思:fastq-dump -X 20000 --split-files -O sra SRR1972877
#* -X 20000是指拆分获取前20000个spot
#* 如果大家了解fastq的格式的话,一个spot就是指一条完整的测序信息,包含四个组成:名字、序列、+号(+号后面可以加别的说明信息)和测序质量信息。
#* --split-files就是把双端测序拆分成两个fastq文件,并且会自动在名字(SRR1972877)后加后缀。
#* -O sra 是指把拆分好的文件放到sra这个文件夹下。
#* 我们可以ls -al查看一下:
bash get-data.sh

# 现在来比对。
# align-ebola 的脚本要先存在。该脚本将比对一个SRA文件并产生一个BAM文件。
#* align-ebola 的脚本在本章节末。

cat sraids.txt | awk ' { print "bash align-ebola.sh " $1 } ' > align-data.sh
bash align-data.sh

# 从所有的数据集中Call snps。
freebayes -p 1 -f ~/refs/852/KM034562.fa bam/*.bam > freebayes.vcf

Alignment script: align-ebola.sh
# 比较两种比对方式的输出。
# 用法: bash align-ebola.sh SRA-RUN-ID
# 把比对文件放到一个单独的文件夹。
mkdir -p bam
# 在错误处停止运行。
set -ue
# 数据文件名。
NAME=$1
# 从传入的参数创建read名字。
# ABC会变为sra/ABC_1.fastq 和 sra/ABC_2.fastq
# 形成输入名称。
READ1=sra/"$NAME"_1.fastq
READ2=sra/"$NAME"_2.fastq
# 参考序列文件。必须两种比对软件都对其进行索引。
REFS=~/refs/852/KM034562.fa
# 我们要加一个read group到mapping
#* \t表示的是tab分隔
GROUP="@RG\tIDNAME\tSMNAME\tLBNAME"
# 运行bwa并生成bwa.sam。
# 生成一个排过序的bam文件
# 这将是bam文件的名字。
BAM=bam/$NAME.bam
bwa mem -R $GROUP $REFS $READ1 $READ2 2> logfile.txt | samtools view -b - | samtools sort -o - tmp > $BAM
echo $?
samtools index $BAM
echo "*** Created alignment file $BAM"
# 打印简单的统计。
samtools flagstat $BAM | grep proper


Lecture 24 - 创建间隔数据

#* 这个lecture,是为后面两个lecture做一些准备
# 获取埃博拉病毒基因组数据和特征(假如你还没有)。
# 一个演示的基因组。也可以是其他的。
efetch -id KM034562 -format gb -db nucleotide > KM034562.gb
readseq -format=FASTA -o ~/refs/852/KM034562.fa KM034562.gb
readseq -format=GFF -o KM034562.gff KM034562.gb

# 我们将会手动用一个编辑器创建BED和GFF文件,并展示它们。
# BED文件示例, 列是tab分隔的:(当你用浏览器复制粘贴,通常会变为了空格)

KM034562    100    200    Happy    0    -

# Example BEDgraph file

KM034562    100    200    50KM034562    150    250    -50


# GTF示例文件
KM034562    lecture    CDS    100    200    .    -    .    gene "Happy"; protein_id "HAP2"

# GFF示例文件
##gff-版本3
KM034562    lecture    CDS    100    200    .    -    .    gene=Happy; protein_id=HAP2

# 数据展示:分割示例
# GTF的
##gff-版本2

KM034562    demo    exon    1000    2000    .    +    .    gene_id "Happy"; transcript_id "Sad";
KM034562    demo    exon    3000    4000    .    +    .    gene_id "Happy"; transcript_id "Sad";
KM034562    demo    exon    5000    6000    .    +    .    gene_id "Happy"; transcript_id "Sad";
KM034562    demo    exon    7000    8000    .    +    .    gene_id "Happy"; transcript_id "Sad";
KM034562    demo    exon    1000    2000    .    +    .    gene_id "Happy"; transcript_id "Mad";
KM034562    demo    exon    5000    6000    .    +    .    gene_id "Happy"; transcript_id "Mad";
KM034562    demo    exon    7000    8000    .    +    .    gene_id "Happy"; transcript_id "Mad";
KM034562    demo    exon    1000    2000    .    +    .    gene_id "Happy"; transcript_id "Bad";
KM034562    demo    exon    5000    6000    .    +    .    gene_id "Happy"; transcript_id "Bad";

# GFF展示
##gff-版本 3

KM034562    demo    exon    1000    2000    .    +    .    Parent=Sad,Mad,Bad;
KM034562    demo    exon    3000    4000    .    +    .    Parent=Sad;
KM034562    demo    exon    5000    6000    .    +    .    Parent=Sad,Mad,Bad;
KM034562    demo    exon    7000    8000    .    +    .    Parent=Sad,Mad;

# BED 展示


KM034562    999    8000    Sad    0    +    999    8000    .    4    1000,1000,1000,1000    0,2000,4000,6000
KM034562    999    8000    Mad    0    +    999    8000    .    3    1000,1000,1000    0,4000,6000
KM034562    999    8000    Bad    0    +    999    8000    .    3    1000,1000    0,4000



欢迎到微信公众号订阅我们
生信媛
bio_sxy




本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:GEO里出现一个GSE里面有多个平台GPL
下一篇:GWAS的困境和遗传模型的新思
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 02:13 , Processed in 0.041744 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.