007-生信人的Linux考试

Jul 2, 2018 10 min read cnposts

刘小泽写于18.7.2 学完了理论，开始实践一下吧主要参考了jimmy提供的练习题

一、在主目录下面创建/tmp文件夹，并且使其中包含 1/2/3/4/5 格式的文件夹系列

mkdir -p  /tmp/1/2/3/4/5
#如果要删除的话，使用rm -r 1/

二、在创建好的文件夹下面，比如我的是 /home/u1239/tmp/1/2/3/4/5 ，里面创建文本文件 xi.txt，并输入内容，例如Hello world, Welcome to bioinfoplanet, Nice to see you。三句话分行显示

echo -e "Hello world,\nWelcome to bioinfoplanet,\nNice to see you" | tee /tmp/1/2/3/4/5/xi.txt 
#利用管道符将echo输出的内容赋给tee命令；tee命令作用是：将接收到的数据一份输出到屏幕，一份到文件

三、在tmp/下创建 1~5这5个文件夹，然后每个文件夹下面继续创建 1~5这5个文件夹, 并查看

mkdir -p /tmp/{1..5}/{1..5} # 然后tree查看

四、增加一点难度：我想在练习三的每个目录中都放进去一个文件xi.txt，内容还是练习二的内容

echo {1..5}/{1..5} |xargs -n1 cp xi.txt  
#删除1~5这5个文件夹及其子文件夹：rm -r {1..5}

五、下载 http://www.biotrainee.com/jmzeng/igv/test.bed 文件，统计该文件总共有几行，含有 H3K4me3 的是第几行

#统计行数：
cat test.bed | wc -l 
#包含H3K4me3: cat test.bed | grep -n "H3K4me3"

六、下载 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件，并且解压，查看里面的文件夹结构

wget  http://www.biotrainee.com/jmzeng/rmDuplicate.zip && unzip *.zip 
#查看结构 tree rmDuplicate

七、打开第六题解压的文件，进入 rmDuplicate/samtools/single 文件夹里面，查看后缀为 .sam 的文件，搞清楚生物信息学里面的SAM/BAM定义是什么【另外关于sam头部注释信息在tmp.header中】

SAM的全称是sequence alignment/map format, 主要用于测序序列mapping到基因组上的结果表示 BAM是SAM的二进制压缩文件，不能像sam可以用less查看，它需要用samtools view打开 SAM分数据通常有两部分，头部注释信息（header section ）和主体比对结果部分（alignment section）
注释信息：SAM文件产生以及被处理过程的一个记录，@开头，包括这么几项：
@HD，说明符合标准的版本、对比序列的排列顺序
@HD VN:1.0 SO:coordinate （排序类型）头部区第一行，VN是格式版本；SO表示比对排序的类型，有unknown（default），unsorted，queryname和coordinate几种。
samtools软件在进行行排序后不能自动更新bam文件的SO值，而picard却可以
@SQ，参考序列说明@SQ
SN:chr1 LN:248956422 （序列ID及长度） SN是参考序列名；LN是参考序列长度；每个参考序列为一行
@RG，比对上的序列（read）说明
@RG ID:sample01 （样品基本信息） 1个sample的测序结果为1个Read Group，每个sample可以有多个文库的测序结果这些内容可以用下面👇程序添加，这些内容可以在形成sam时加入，其中ID是必须要有的 bwa mem -R “@RG\tID:$sample\tSM:$sample\tLB:WES\tPU:Illumina\t PL:Miseq" 【ID：样品ID号；SM：样本名；LB：文库名；PU：测序仪；PL：测序平台】【注意一定要在@RG后加入\t，字段之间用tag分割，否则运行错误】
@PG，使用的程序说明 @PG ID:bowtie2 PN:bowtie2 VN:2.2.9
@CO，任意的说明信息
主体部分：11个主列和1个可选列，之间用tab分割【如果某一列为*或者0，说明这一列没有信息】
第一列：QNAME—比对的序列名称 例如：SRR1042600.42157053 ，也就是read ID 如果reads比对到多条序列或者某条序列的多个位置，相同的名字会出现多次。【有可能read比对到两个位置，mate也比对到两个位置，这样这一条reads就有四次相同名字出现】 Illumina有两种双端测序方式：PE测序和MP测序。 MP测序中：read表示本条read，mate表示pair中的另一条read；或者在PE测序中用read1，read2表示
第二列：Flag read mapping情况的数字表示，较少的数字记录表示当前序列的匹配情况 1: 表示序列PE双端测序，read是序列的一条 2: read和参考序列完全比对，没有插入缺失等 4: 这个序列没有比对上 8: 这个序列的另一端没有比对上 16: 这个序列比对到参考序列负链上 32: 这个序列的另一端比对到参考序列负链上 64: 表示这个序列是read / read1 128: 表示这个序列是mate / read2 256: 这个序列是次优的比对结果，它是secondary 【一条序列可能比对到参考序列的多个位置，只有一个是首要的比对位置】 512:序列QC时失败，被过滤掉 –> 不常用 1024：PCR或测序错误产生的重复reads –> 不常用 2048: 补充匹配的序列 –> 不常用
比如flag 99=64+32+2+1, 意思就是：来自于双端测序的read1，且匹配的非常好，对应的read2匹配到了负链(该read1为正链)；
83=64+16+2+1表示read1完全匹配负链，
而147=128+16+2+1表示双端测序的read2完全匹配负链
163=128+32+2+1 表示read1完全匹配负链
161=163-2,81=83-2. 表明这些read不是完全匹配，存在插入缺失
        
#常见的FLAGs:
其中一条reads没有map上: 73, 133, 89 121, 165, 181, 101, 117, 153, 185, 59, 137
【PE测序的两条read可能同时出现上面的数字，可能表示一条没有匹配上，另一条完美匹配，也可能都不匹配。因此这个数字并不能精确告诉你，哪条read是匹配上的，只能模糊的说，他判断出有这种情况，要再结合CIGAR判断】
两条reads都没有map上: 77,141
比对上了，方向也对，也在插入大小(insert size)内: 99, 147, 83, 163
比对上了，也在插入大小(insert size)内， 但是反向不对:67, 131, 115, 179
单一配对，就是插入大小(insert size)不对: 81, 161, 97, 145, 65, 129, 113, 177
第三列：Rname read比对的那条序列的序列名称，与头部注释文件中@SQ中的SN对应；没有比对上用*
第四列：POS read比对到RNAME这条序列的最左边的位置（从1开始计数，没有比对上，此处为0）
第五列：MAPQ mapping的质量值，即这段序列匹配到参考序列某段的概率值，结果= -10log10(Err) , 一般为0-60，60表示比对率最高假如一条read比对到了两处，但是Q值不同，一般Q值大的那个比对的真实性更高
第六列：GIGAR reads比对到Rname的具体情况，使用数字+字母表示
CIGAR 值 BAM 含义
M 0 alignment match (can be a sequence match or mismatch)
I 1 insertion to the reference
D 2 deletion from the reference
N 3 skipped region from the reference
S 4 soft clipping (clipped sequences present in SEQ)
H 5 hard clipping (clipped sequences NOT present in SEQ)
P 6 padding (silent deletion from padded reference)
= 7 sequence match
X 8 sequence mismatch
S和H一般用于read前后出现大部分的错配，但是中间能够比对上的情况 H只出现在CIGAR的首末尾；S与H在CIGAR中一般组合出现，但也有只存在S的情况 S表示序列会出现在SEQ中，H则不会出现在SEQ列中 M与=的差别：M能配对，但可能存在某个碱基的差异，=表示完全匹配对于mRNA比对到基因组这种情况，N的出现代表内含子存在 M/I/S/=/X这5项的加和应该与序列长度相等
【以下第7-9列是mate的信息，如果是单端测序，这几列没有意义】 第七列：RNEXT 双端测序中read2比对到的参考序列名称【也就是mate比对到的染色体号】=表示read1、read2都比对到了同一条；*表示没有序列名称；如果read2匹配到的和read1不同，那么RNEXT会更换名称
第八列：PNEXT
read2/mate比对到参考序列上的第一个碱基位置【read1的是第四列的POS】
第九列：TLEN 估计的模板片段长度，就是read寻找到了参考序列的哪些区域；当mate/read2位于本序列上游时(就是位于read1之前)，该值为负值； PE测序中，利用这个TLEN值可以推测出文库的平均大小
第十列：segment序列
第十一列：segment质量值
可选列： 格式一般为TAG:TYPE:VALUE ,例如：NM:i:4
NM:i 经过编辑的序列长度(edit distance)
MD:Z 错配位置/碱基(mismatching positions/bases)
AS:i 匹配得分(Alignment score)
XS:i 第二好的匹配得分(suboptimal alignment score)
YS:i mate序列匹配的得分
BC: 条码序列（barcode sequence)

CIGAR 值	BAM	含义
M	0	alignment match (can be a sequence match or mismatch)
I	1	insertion to the reference
D	2	deletion from the reference
N	3	skipped region from the reference
S	4	soft clipping (clipped sequences present in SEQ)
H	5	hard clipping (clipped sequences NOT present in SEQ)
P	6	padding (silent deletion from padded reference)
=	7	sequence match
X	8	sequence mismatch

八、bam文件的生成与查看？

#使用bowtie2进行mapping，bowtie2速度上比较有优势，而且输出的结果就是SAM格式的
#1.建立fa文件索引,以人类基因组hg19为例,设定前缀为hg19
bowtie2-build hg19.fa hg19 --> 将会生成好几个前缀hg19，后缀bt2的文件
#2.将PE测序数据比对到基因组
bowtie2 -q --phred64 -p8 --no-unal -x hg19 -1 read1.fq.gz -2 read2.fq.gz -S paired.sam
#-q:输入格式为fq
#--phred 33设定碱基质量值（默认33，现在一般测得都是33）
#-p:线程数
#--no-unal:将没有比对上的剔除(--no-unalignment)
#-x：bowtie2的index参考序列
#-S：设置输出sam名称

#3. 使用samtools将sam转bam
samtools view -b -S file.sam > file.bam
#还能将bam转sam
samtools view -h file.bam > file.sam

#4. bam文件读取查看
samtools view xxx.bam |less 
#如果只想查看header信息，就加一个-H参数 (-h:显示header信息)

九、查看使用的参考基因组有多少染色体

第八题中，我们使用了hg19基因组，用bowtie2比对得到sam，又转为bam
因此，我们直接对bam进行操作
操作流程：查看bam的头注释文件中的第二列SN:chrxxx，其中包含了染色体信息--> sed去掉数字前的信息 -->对数字排序
samtools view -H tmp.rmdup.bam | cut -f2 | sed -e 's/SN:chr//g' | sort -n

十、rmDuplicate/samtools/single中的sorted.bam文件中，统计第二列各flag个数

samtools view tmp.sorted.bam | cut -f2 | sort -n | uniq -c

同理也可以在rmDuplicate/samtools/paired中的sorted.bam文件中统计flag

**十一、下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件，并且解压，进入解压后的文件夹test1_fastqc，找到 fastqc_data.txt 文件，统计以»开头的有多少行 **

cat fastqc_data.txt | grep '>>' | wc -l
# 结果会显示24行

十二、下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss (存储转录起始位点信息)文件，在NCBI中找到自己感兴趣的基因对应的 refseq数据库 ID（这里我以乳腺癌BRCA1为例），然后找到它在hg38.tss 文件的哪一行

#这里搜索的是第一个RefSeq ID
cat hg38.tss | grep 'NM_007294' 
#得到结果 NM_007294	chr17	43123483	43127483	1

关于RefSeq数据库：详细的会在之后的更新中介绍，关于各个数据库的使用豆豆和花花会整理
RefSeq数据库中所有的数据是一个非冗余的、提供参考标准的数据，包括染色体、基因组（细胞器、病毒、质粒）、蛋白、RNA等。refseq序列是NCBI筛选过的非冗余数据库，比GeneBank准确。 **关于它的ID：**NM开头的表示标准序列，XM表示预测的蛋白编码序列，NR表示非编码蛋白的mRNA序列，AF开头的表示克隆序列，BC开头的表示模板序列另外，你可能见过gi|4557284|ref|NM_000646.1|[4557284]这种格式 gi就是代表genebank identifier；ref就是对应的refseq中的ID啦

十三、解析hg38.tss 文件，统计每条染色体的基因个数

cut -f2 hg38.tss |cut -d’_’ -f1 | sort |uniq -c | sort -rn

 6157 chr1
   5880 chr19
   5782 chr6
   4090 chr2
   3794 chr17
   3577 chr11
   3395 chr3
   3014 chr12
   2838 chr10
   2821 chr5
   2785 chr7
   2696 chr16
   2561 chrX
   2377 chr15
   2310 chr9
   2277 chr4
   2221 chr8
   1982 chr14
   1692 chr20
   1410 chr22
   1133 chr13
    895 chr21
    883 chr18
    414 chrY
     32 chrUn
      2 chrM
#chrM是线粒体染色体
#chrUn是不能准确地放到特定染色体上的基因组DNA克隆重叠群

十四、统计hg38.tss 文件中NM和NR开头的序列，了解NM和NR开头的含义

#统计NM开头
cat hg38.tss | grep 'NM' | wc -l 
#同理可以统计NR开头
NM开头的表示标准序列，可以转录成蛋白质的基因
NR非编码蛋白的mRNA序列

linux practice

007-生信人的Linux考试

Yunze Liu

Bioinformatics Sharer

Related