Subscribed unsubscribe Subscribe Subscribe

CPSAM.org

computer, programming, statistics and more!

linuxファイルサーバ

sambaをインストール /var/sambaを\hostaname\data として認識させるためのsmb.confファイルの例

smb.conf

[global] workgroup = WORKGROUP server string = %h server (Samba, Ubuntu) dns proxy = no log file = /var/log/samba/log.%m max log size = 1000 syslog = 0 panic action = /usr/share/samba/panic-action %d server role = standalone server passdb backend = tdbsam obey pam restrictions = yes unix password sync = yes passwd program = /usr/bin/passwd %u passwd chat = Enter\snew\s\spassword: %n\n Retype\snew\s\spassword: %n\n password\supdated\ssuccessfully . pam password change = yes map to guest = bad user

[data] path = /var/samba writeable = yes read only = No force create mode = 0666 force directory mode = 0777

Remove all ads

RNA-seq pipeline

リファレンスゲノムをダウンロードする
併せてアノテーションファイルもダウンロードする

wget ftp://ftp.ensembl.org/pub/release-85/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz
wget ftp://ftp.ensembl.org/pub/release-85/gtf/homo_sapiens/Homo_sapiens.GRCh38.85.gtf.gz

bowtieでindexing

bowtie2-build -r Homo_sapiens.GRCh38.dna.toplevel.fa \
  Homo_sapiens.GRCh38.dna.toplevel

Tophatでgtfファイルもindexingしておく
下記の例だとカレントディレクトリにtranscriptomeというディレクトリが作成され
その中にHomo_sapiens.GRCh38.85_gtfというプレフィクスでいくつかのファイルが作成される

tophat -G Homo_sapiens.GRCh38.85.gtf \
  transcriptome/Homo_sapiens.GRCh38.85_gtf \
  Homo_sapiens.GRCh38.dna.toplevel

サンプルデータをダウンロードする

sra=SRR213838
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP061/SRP061888/${sra}{1..8}/${sra}{1..8}.sra

sratoolkitを使ってsraファイルを解凍する。 --split-3オプションを使わないと、ペアエンドデータに解凍できないので注意。

fastq-dump --split-3 ${sra}{1..8}.sra

FastQCにかける

fastqc --nogroup ${sra}{1..8}_1.fastq
fastqc --nogroup ${sra}{1..8}_2.fastq

アラインメント

tophat -p 2 \
       -G Homo_sapiens.GRCh38.85.gtf \
       Homo_sapiens.GRCh38.dna.toplevel \
       SRR2138381_1.fastq \
       SRR2138381_2.fastq
Remove all ads

fasta整形

簡単にいうとbiomartからとってきたutr sequenceが改行されていて大変だったという話

改行の置換はsedは苦手らしい

cat martquery.txt | tr -d "\n" > out

sed -E "s/>/\n>/g" out >  out2
sed -E "s/([0-9]){1}([AGCTS])/\1,\2/g" out2 > out3
sed -E "s/^>//" out3 > out4
sed -E "s/\|/,/" out4 > utr.csv
Remove all ads

Box同期フォルダの変更

box-syncのインストール前にregeditで HKEY_LOCAL_MACHINE/SOFTWARE/Box/BoxSync にString variableを追加、名前はSyncRootFolder、値をD:\boxなどとして保存 box-syncを起動、ログインで指定したディレクトリが作成 ※このときにすでにdirectoryがあるとエラーになる

Remove all ads

マイナスから始めるExome解析

0から始めるエクソームデータ解析
を参考にさせていただきました

とりあえず現行では

sudo apt-get install ant
sudo apt-get install gradle
git clone https://github.com/broadinstitute/picard.git
cd picard
git clone https://github.com/samtools/htsjdk.git
./gradle -jar
cd ../
ant

でpicardのbuildまでは成功する


Broad instituteからbudle(reference genome, interval fileなどがセットになったもの)の入手
ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b2.8/b37より適宜入手する

Reference genome: human_g1k_v37.fasta, human_g1k_v37.fasta.fai
Variant file: dbsnp_138.b37.vcf, dbsnp_138.b37.vcf.idx
Insert-deletion file: 1000G_phase1.indels.b37.vcf, 1000G_phase1.indels.b37.vcf.idx


サンプルデータの入手(http://www.ncbi.nlm.nih.gov/pubmed/21936905)

wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/SRA044/SRA044780/SRX091464/SRR330441_1.fastq.bz2
wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/SRA044/SRA044780/SRX091464/SRR330441_2.fastq.bz2
bunzip2 *.fastq.bz2

Reference genomeにindexをつける

bwa index -p human_g1k_v37.fasta -a bwtsw human_g1k_v37.fasta

picardでreferenceのdictionaryを作る

java -jar ../tools/picard.jar CreateSequenceDictionary \
 R=human_g1k_v37.fasta \
 O=human_g1k_v37.dict

アラインメント

bwa aln -t 2 human_g1k_v37.fasta SRR330441_1.fastq > SRR330441_1.sai
bwa aln -t 2 human_g1k_v37.fasta SRR330441_2.fastq > SRR330441_2.sai
bwa sampe -P human_g1k_v37.fasta \
 -r '@RG\tID:01\tSM:s6\tPL:Illumina' \ 
 -f SRR330441.sam \
 SRR330441_1.sai SRR330441_2.sai \
 SRR330441_1.fastq SRR330441_2.fastq

bamへの変換とソート、インデックス付与

samtools view -bS SRR330441.sam > SRR330441.bam
samtools sort SRR330441.bam -o SRR330441.sorted.bam
samtools index SRR330441.sorted.bam

Picardによる重複リードの除去とこれに対するインデックス付与

java -Xmx2G -jar ../tools/picard.jar MarkDuplicates \
 ASSUME_SORTED=true \
 REMOVE_DUPLICATES=true \
 INPUT=SRR330441.sorted.bam \
 OUTPUT=SRR330441.removed.sorted.bam \
 METRICS_FILE=duplicate \
 VALIDATION_STRINGENCY=LENIENT

samtools index SRR330441.removed.sorted.bam


bamファイル名を記載した"files.list"を作る

echo 'SRR330441.removed.sorted.bam' > files.list

GATKのRealignerTargetCreatorを使って既知のinsert/deletion付近の再マッピング用のファイルを作る

java -jar ../tools/GenomeAnalysisTK.jar \
 -T RealignerTargetCreator \
 -R human_g1k_v37.fasta \
 -I files.list \
 --known ../bundle/1000G_phase1.indels.b37.vcf \
 -log intervals.log \
 -L ../sureselect/S07604514_AllTracks.bed
 -o SRR330441.intervals
Remove all ads

Ubuntu server 16.04

Ubuntu Server 16.04にシステムを入れ替えたところ画面が真っ暗

いろいろ調べると、recoveryモードで立ち上げて

>||

vi /etc/default/grub

GRUB_CMDLINE_LINUX=""

-> GRUB_CMDLINE_LINUX="vga=773"

GRUB_CMDLINE_LINUX_DEFAULT=""

-> GRUB_CMDLINE_LINUX_DEFAULT="text"

update-grub

shutdown -r now

||<

これで解決

Remove all ads

sed

>表示

sed -n 1p #一行目を表示
sed -n 1~2p #1+2n行目、つまり奇数行を表示
sed -n 1~5p #1+5n行目、つまり1,6,11,16...行目を表示
sed -n -E /foo/p #正規表現にマッチする行を表示

>置換
sed -E "s/foo/bar file.txt" #最初に見つかったfooをbarに置換
sed -E "s/foo/bar/3 file,txt" #3番目に見つかったfooをbarに置換
sed -E "s/foo/bar/g file.txt" #行内のfooをbarに置換

>検索
sed -n -E "/foo/p file.txt" #fooの含まれる行のみ表示
sed -n -E "/foo/d file.txt" #fooの含まれない行のみ表示

正規表現を使うときは””で囲む

Remove all ads